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^ \ The most promising source of gravitational waves for the planned kilometer-size 

laser-interferometer detectors LIGO and VIRGO are merging compact binaries, i.e., 
neutron star/neutron star (NS/NS), neutron star /black hole (NS/BH), and black 
hole/black- hole (BH/BH) binaries. We investigate how accurately the distance to 

o 

g : the source and the masses and spins of the two bodies will be measured from the in- 

\ spiral gravitational wave signals by the three detector LIGO/VIRGO network using 

qh| "advanced detectors" (those present a few years after initial operation). The large 

^j-jl number of cycles in the observable waveform increases our sensitivity to those param- 



> 



eters that affect the inspiral rate, and thereby the evolution of the waveform's phase. 



H ' These parameters are thus measured much more accurately than parameters which 

affect the waveform's polarization or amplitude. To lowest order in a post-Newtonian 
expansion, the evolution of the waveform's phase depends only on the combination 
M = (MiM2)^/^(Mi + M2)~^/^ of the masses Mi and M2 of the two bodies, which 
is known as the "chirp mass." To post-l-Newtonian order, the waveform's phase also 
depends sensitively on the binary's reduced mass /i = M1M2/ {Mi +M2), allowing, in 
principle, a measurement of both Mi and M2 with high accuracy. We show that the 
principal obstruction to measuring Mi and M2 is the post-1. 5-Newtonian effect of the 
bodies' spins on the waveform's phase, which can mimic the effects that allow ^ to be 
determined. The chirp mass is measurable with an accuracy /S.^A/M ~ 0.1% — 1%. 
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Although this is a remarkably small error bar, it is ~ 10 times larger than previous 
estimates of ^M/M which neglected post-Newtonian effects. The reduced mass 
is measurable to ~ 10% - 15% for NS/NS and NS/BH binaries, and - 50% for 
BH/BH binaries (assuming IOM0 BH's). Measurements of the masses and spins are 
strongly correlated; there is a combination of /x and the spin angular momenta that 
is measured to within ~ 1%. Moreover, if both spins were somehow known to be 
small (< O.OlMf and < O-OlMf , respectively), then jj, could be determined to within 
~ 1%. Finally, building on earlier work of Markovic, we derive an approximate, 
analytic expression for the accuracy AD of measurements of the distance D to the 
binary, for an arbitrary network of detectors. This expression is accurate to linear 
order in l/p, where p is the signal-to-noise ratio. We also show that, contrary to 
previous expectations, contributions to AD/D that are nonlinear in 1/p are signifi- 
cant, and we develop an approximation scheme for including the dominant of these 
non-linear effects. Using a Monte-Carlo simulation, we estimate that distance mea- 
surement accuracies will be < 15% for ~ 8% of the detected signals, and < 30% for 
~ 60% of the signals, for the LIGO/VIRGO 3-detector network. 
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I. INTRODUCTION 



Neutron star-neutron star (NS-NS) binaries with orbital periods of less than half a day 
will spiral together and merge in less than a Hubble time, due to gravitational radiation 
reaction. Three such short-period NS-NS binaries have been observed in our Galaxy; when 
extrapolated to the rest of the Universe these observations result in an estimated NS-NS 
merger rate in the Universe of ~ 10^ yr^^Gpc^^ 0,0]. A strong gravitational wave signal 
is emitted during the last few minutes of inspiral, before the tidal-disruption/coalescence 
stage begins. If the Laser Interferometer Gravitational Wave Observatory (LIGO) [Q, and its 
French-Italian counterpart VIRGO @], achieve the so-called "advanced detector" sensitivity 
level of Ref. P], then they will be able to detect gravitational waves from the last few minutes 
of NS-NS inspirals out to distances of order ~ 1 Gpc 0. Hence, event rates of order 10^ yr~^ 
may be achieved. While there is no direct observational evidence relevant to the merger rates 
for neutron star-black hole (NS-BH) and black hole-black hole (BH-BH) binaries, arguments 
based on progenitor evolution scenarios suggest that these merger rates may also be on the 
order of 10^ yr-^Gpc-^ [|,|. The merger of two IOMq black holes would be detectable by 
LIGO/ VIRGO out to cosmological distances at redshifts of ~ 2 — 3. 

The gravitational waveforms arriving at the detectors depend on the inspiraling bodies' 
masses and spins, the distance to the binary, its angular position on the sky, and the orienta- 
tion of the binary's orbital plane. By comparing the observed waveforms with theoretically 
derived templates, the observers will extract these parameters to a level of accuracy that 
is determined by the noise in the detectors, and by the detectors' relative positions and 
orientations. From the output of a single detector, there will be sufficient information to 
determine the masses of the two bodies, but not their distance or their location on the sky. 
By combining the outputs of the three LIGO/ VIRGO detectors, it should be possible to de- 
termine the location of the binary on the sky to within ~one degree PJ^, and the distance 
to the binary to within ~ 30%. 

There are many potential applications of such measurements, as has been emphasized 
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by Schutz P]. For example, coalescing binaries are potentially very useful standard candles 
for astronomical distance measurements — it has been estimated that from ~ 10^ detected 
NS-NS events it will be possible to determine the Hubble constant Hq to within ~ 10% 
PJ9| jTO|JTl[] . It may also be possible to measure NS radii, and thus constrain the NS equation 
of state, by measuring the frequency at which the tidal disruption of the neutron star causes 
the waves to shut off [§]. And from gravitational wave observations of the final coalescence 
of two black holes, there may follow new insights into gravitational dynamics in the highly 
nonlinear regime. The effectiveness of these and other applications depends on the accuracy 
with which one can read off, from the measured waveform, parameters such as the distance 
to the binary and the masses of its two components. 

The purpose of this paper is to estimate the limits on measurement accuracies that arise 
from sources of noise that are intrinsic to the detectors. These sources of noise include, for 
example, thermal vibrations in the interferometers' suspended test masses, and randomness 
in the arrival times of individual photons at the interferometers' mirrors (photon shot- 
noise), which simulate in the interferometers' output the effects of gravitational waves 0. 
Intrinsic detector noise is expected to be the dominant source of error in the determination of 
coalescing binary parameters, in part because gravitational waves interact very weakly with 
matter through which they pass ||12|. Other possible sources of error which we do not consider 



here include (i) systematic errors due to insufficiently accurate theoretical modeling of the 
gravitational waveforms, which will be important primarily for mass and spin measurements 



||8|JT3[|; and (ii) amplification/deamplification of the wave amplitudes by gravitational lensing 
effects, which will be important primarily for distance measurements. See Markovic for 
a detailed discussion of this issue. 

Many of our conclusions have already been summarized in Cutler et al. M. Initial 



measurement accuracy analyses have been carried out by Finn and Chernoff |T^, and by 
Jaranowski and Krolak |[T^ , using a simplified, "Newtonian" model of the waveform. While 
Newtonian waveforms are adequate for predicting how accurately one can measure the dis- 
tance to the source, they do not allow one to calculate how accurately the individual masses 



can be measured, as we explain below. For this purpose, one must include post-Newtonian 
corrections to the waveform. 

Much of our work was guided by the following sequence of considerations. These con- 
siderations introduce some of the issues addressed in this paper, motivate a number of 
approximations that we make in our analysis, and give a preview of some of our main 
conclusions. 

First, coalescing binaries are very "clean" sources of gravitational waves: the waveform 
is determined to high accuracy by a relatively small number of parameters These pa- 
rameters are the source's location, orientation, time of coalescence, and orbital phase at 
coalescence, as well as the bodies' masses and spin angular momenta. Various other compli- 
cating physical effects, not described by these parameters, can be shown to be unimportant. 
We can generally assume, for instance, that the orbits are circular. This is because radia- 
tion reaction causes the orbit's eccentricity e to decrease during the inspiral, according to 
e"^ oc pi9/9, where P is the orbital period [|16|. (The effect of a small eccentricity on the 



phase of the waveform scales like e^.) The emitted gravitational waves are in the frequency 
band accessible to LIGO only for the last few minutes of inspiral, when P < 0.2 sec. Thus 
a binary born with e of order unity and P > 1 hour will have < 10^^ by the time it be- 
comes "visible" to LIGO ll^. Also, tidal interactions between the bodies have been shown 
to be negligible |ll8| , p!9| (except for the last few orbits), so for our purposes the bodies can 



be treated as structureless, spinning point masses pO| 



Second, this high predictability of the gravitational waveforms means that the technique 



of matched filtering can be used to detect the waves [|T2|. For the most distant (most 
frequently observed) sources, this will involve extracting the waveforms from the considerably 
larger instrumental noise in which they will be imbedded. The technique works as follows 



r^ . The measured strain amplitude in each detector 

sit) = hit) + n{t) (1) 
consists of a (possibly present) signal h{t), and the detector noise n(t), which we assume is 



Gaussian. To detect any imbedded signal, one first suppresses those frequency components 
of tlie signal at which the detector noise is largest by convolving with Wiener's optimal filter 
wit): thus, s{t) — > / w{t — r)s{T)dT |21|]. Then, for each inspiral waveform h(t) in a large 
set of theoretical template waveforms, one computes the signal-to-noise ratio S/N, defined 
by 

S_^^ ^ IHt)w{t-T)s{T) drdt 
N rms / h{t) w{t — r) n{T) drdt 

In Eq. (^, the denominator is what would be the root-mean-square value of the numerator, 
if the detector output ([l|) consisted of noise alone. Thus, when no gravitational wave is 
present, each S/N[h] is a random variable with Gaussian distribution and root-mean-square 
equal to 1. Conversely, if S/N[h] is sufficiently large as to basically preclude the possibility of 
its arising from noise alone — for any of the ~ 10^^ template waveforms that will be applied 
to the data each year — then one can assert with high confidence that a gravitational wave 
h has been detected, and that h is close to h. It is easy to show that if some template 
waveform h yields a signal-to-noise ratio of S/N[h] > 6.0 in each of two detectors, then 
with > 99% confidence a gravitational wave has been detected ||^. Defining the combined 
signal-to-noise ratio p of a network of detectors by 



P^^T^Pl (3) 

where pa is the S/N in the ath detector, we see that p ~ 8.5 represents the "detection 
threshold" for two detectors. For a three-detector network, the detection threshold is still 
p ~ 8.5, corresponding to S/N ^ 4.9 in each detector. Since detections at threshold 
represent the most distant coalescences that one can observe (given the binary's masses, its 
orientation, and its angular position on the sky), and since coalescing binaries are presumably 
distributed roughly uniformly on large scales 100 Mpc), the mean value of p for detected 
events will be roughly 1.5 times the threshold value |2^. Thus "typical" detections will have 
p ^ 12.7. Similarly, the strongest 1% of signals should have p ^ 40; i.e., (100)^/^ times the 
threshold value. 
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Third, much more information is obtainable from the waveform than one might naively 
expect, for the following reason. The LIGO and VIRGO detectors will be broad-band 
detectors, with good sensitivity in the frequency range 10 — 500 Hz. The gravitational wave 
trains from inspiraling stellar-mass binaries typically contain ~ 10'^ cycles in this range. 
Now, if the signal h{t) and template h{t) lose phase with each other by just one cycle out of 
thousands, as they sweep upwards in frequency from ~ 10 Hz to ~ 500 Hz, then the integral 
/ h{T)h{t)w{t — T)dTdt will be significantly diminished. Consequently the value of S/N[h] 
will be small unless the phase of the template waveform h is "just right" throughout the 
inspiral. Since the evolution of the waveform's phase is largely determined by the masses of 
the two bodies (through their influence on the inspiral rate), one might expect to measure 
the masses of the bodies with fractional error ~ l/A/'cyc, where A/'cyc is the total number of 
cycles in the observed waveform. This fractional error of ~ 10^'^ contrasts with the ~ 20% 
accuracy with which one can determine parameters, such as the distance to the source, that 
do not affect the phase evolution (as was first pointed out by Cutler et al. and by Chernoff 
and Finn [|T1|). 



Fourth, our extension of the measurement-error analysis to include post-Newtonian ef- 
fects introduces the following new features. To Newtonian order, the gravitational wave 
signal depends on the two masses only through the particular combination Ai = fx^^^M"^^^, 
where /i is the reduced mass and M is the total mass of the system. This combination is 
referred to as the "chirp mass." The degeneracy in the dependence on the masses is broken, 
however, by post-Newtonian effects that in principle allow one to determine the individ- 
ual masses Mi and M2. In the equation governing the evolution of the waveform's phase 
[Eq. (^) below], the post-Newtonian terms are ~ M/r times smaller than the Newtonian 
terms, where r is the orbital separation. Since M/r ^ 1/20 when the signal is strongest, one 
might expect to determine each of the two masses ~ 20 times less accurately than Ai. We 
show in Sec. |IT1| that this expectation is correct, provided the spins of the bodies are known 
to be small. 

Now, black holes and neutron stars in merging binaries may or may not be rapidly 



spinning. However if we cannot assume a priori that their spin angular momenta are very 
small, then in attempting to find the best fit to the data, we must allow for the possibility that 
the spin angular momenta are of order their maximum possible values. We show in Sec. [111B| 
that the extra "confusion" introduced by the spin-dependence of the waveform worsens the 
accuracy of individual mass measurements by more than an order of magnitude. This is 
easy to understand: the leading order spin terms in the orbital evolution equation [Eq. (^91) 
below] are only one-half post-Newtonian order higher than the leading terms responsible 
for splitting the mass degeneracy. Therefore the effect on the gravitational waveform of 
errors in Mi and M2 that keep Ai fixed can be approximately masked by somewhat larger, 
compensating errors in its spins. Hence the measured values of masses and spins will have 
strongly correlated errors [cf. Fig. ^ below], thereby increasing mass-measurement errors 
pO| . Our results for measurement accuracies are summarized in Tables || and || and Fig. ^ 
below. 

The rest of the paper is organized as follows. In Sec. ^ we review the anticipated detec- 
tor noise levels, the basic elements of signal processing, and the lowest-order, "Newtonian" 
waveforms. In Sec. |T| we calculate expected mass-measurement accuracies, taking post- 
Newtonian effects into account. We do this in two stages: first neglecting spin effects in 



Sec. |I1I A| , then including them in Sec. [HI B| . Our emphasis is on learning roughly what 



accuracies can be expected — in part because to treat the parameter-estimation problem in 
full generality would be extremely complicated. Therefore, we focus on a somewhat simpli- 
fied "model" of the gravitational waveform, which nevertheless incorporates the effects that 
are most important for determining the mass-extraction accuracy. A further approximation 
which we make is to use a linear error-estimation formalism, which is valid when the errors 
are small (or equivalently, when the signal-to- noise ratio is large). 

Most of the information that allows one to measure the binary masses is contained in 
the phase evolution of the waveform (rather than in its amplitude or polarization). Since 
all detectors in a detector network measure very nearly the same phase evolution, for the 
purpose of estimating mass measurement accuracies, to a good approximation it is adequate 
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to model measurements made by single detector. The mass measurement errors for 
detectors are roughly those for a single detector, divided by -\/iV. 

When measuring the distance D to the binary, on the other hand, one must also de- 
termine the position of the source on the sky and the amplitude and polarization of the 
waveform. Hence, to estimate distance measurement accuracies, we must model measure- 
ments by an entire detector network. However in this case it is a needless complication to 



use post-Newtonian waveforms; as we show in Sec. ^\ and Appendix C below, to a good ap- 
proximation it is adequate to use Newtonian waveforms in the analysis. This is our approach 
in Sec. |V|, where we estimate the distance measurement accuracy /\D attainable by an ar- 
bitrary network of detectors. Jaranowski and Krolak have numerically calculated in several 
specific cases the distance measurement accuracy one can achieve with the LIGO/VIRGO 



network [|^. We provide a greatly simplified, analytic solution to the distance-accuracy es- 



timation problem, using an approximation due to Markovic |^ . The approximation consists 
in neglecting the effect on distance measurement errors of the relatively small uncertainty 
in the angular position of the source on the sky. We derive a relatively simple formula for 
the rms distance error AD in this approximation, which applies to any number of detectors 
with arbitrary orientations. 

This formula is derived using the linear error-estimation formalism mentioned above, and 
consequently is accurate only to linear order in 1/D. We show that, contrary to previous 
expectations, effects which are nonlinear in 1/D have a significant effect (i.e., factors ^ 2) 
on the predicted distance-measurement accuracies, and develop an approximate method of 
calculation which gives rough estimates of these nonlinear effects. This method is based on a 
Bayesian derivation of the (non-Gaussian) probability distribution for the distance D, which 
incorporates our a priori knowledge as well as the information obtained from a gravitational 
wave measurement. The method also allows us to estimate values of AD for binaries that 
are seen nearly face-on, for which, as pointed out by Markovic 0, the linear error-estimation 
method breaks down. Our results for nearly face-on binaries are typically factors of order 2 
to 3 smaller than the upper-limit estimates given by Markovic 0. 



In Appendix A we extend the treatment of signal processing given in Sec. || to incorporate 
(i) an arbitrary number of detectors, (ii) the effects of a |)non knowledge, and (iii) estimation 
of measurement errors beyond the linear, Gaussian approximation. These extensions are 
required in Sees. [llllj| and |V|. We also develop other tools which should be useful in 
future analyses of LIGO/VIRGO measurement accuracies: we derive an expression for the 
minimum signal-to-noise ratio {S/N)^\^ necessary in order that the Gaussian approximation 
for estimation of measurement accuracy be valid, and explain how to treat degenerate points 
in parameter space at which the Gaussian approximation breaks down. 

In this paper we will focus on three fiducial types of binary — NS-NS, BH-NS, and BH- 
BH — with fiducial masses Mbh = IOMq and M^s = 1-4M0 (unless otherwise specified). 
Throughout we use units where G = c = 1. Thus all quantities are measured in units of 
seconds, except where, for convenience, we use units of solar masses. The conversion factor 
is IMq = 4.926 X lO'Sec. 

II. DETECTION AND MEASUREMENT 
OF GRAVITATIONAL WAVEFORMS 

A. Detector characteristics 

In order to decide what information can be extracted from gravitational waveforms, one 
must have a realistic model of the detector noise n{t). This noise will have both Gaussian and 
non-Gaussian components. We will restrict our analysis to statistical errors due to Gaussian 
noise. It is likely that the effects of the non-Gaussian components will be unimportant due 
to (i) the rejection of events that are not simultaneously detected in two or more detectors, 
and (ii) the filtering of the detector outputs with theoretical waveform templates; however 
this issue needs further study. 

The remaining Gaussian noise can be described by its spectral density Sn{f), where / is 
frequency. The LIGO team has published an estimate of the noise spectrum that might be 
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attained a few years after LIGO comes on line — the so-called "advanced detector" noise 
spectral density We use the following rough analytic fit to their noise curve: 

loo /<10Hz, 

^n(/)= (4) 

[5o[(/o//r + 2(l + (/V/2))] />10Hz 

where 5*0 = 3 x 10"''^Hz~^ and /o = 70 Hz. For frequencies / < 10 Hz, the noise due 
to seismic vibrations is so large that we take it to be effectively infinite. Thermal noise 
dominates in the frequency band lOHz ^ / ^ 50Hz, and photon shot noise dominates for 
/ ^ 50 Hz. We refer the reader to Refs. |p|,[l2| for more details on the sources of noise. 

The amount of detector noise determines the strength of the weakest signals that can 
be detected, and thus the distance to which a given type of source can be seen. The noise 
level (^) will permit the detection of NS-NS mergers out to ~ 1 Gpc 0,14,^], giving an 



estimated detection rate of ~ 10^ yr^^ 0. In this paper we are principally concerned not 
with detection issues, but rather with the accuracy of parameter estimation. This accuracy 
of parameter estimation depends only on the shape of the noise spectrum, and on the signal- 
to- noise (S/N) of the detection; i.e., simultaneously doubling both the noise levels and the 
signal strength leaves measurement accuracy unchanged. We normalize our results to a 
fixed S/N, and hence our results are independent of the parameter 5*0 appearing in the 
noise spectrum (^. 

Since the LIGO team's publication P| of their estimate of the advanced detector's noise 
curve, there have been new developments in the understanding of the detector's thermal noise 
which indicate that the advanced thermal noise spectrum may be flatter than previously 
thought [^. A modified noise-curve estimate, reflecting this new understanding, has not 
yet been published. Like the noise spectrum in Ref. on which our simplified model (^) is 
based, the modified noise curve will depend on the values of advanced detector parameters 
(such as the quality-factors of modes of vibration of the suspension wires and suspended 
masses) for which only rough estimates are available. Our approximate analytic formula 
describing the modified advanced detector noise curve, 
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Snif) 



OO /<10Hz, 
Sma-^f/U)-' 10Hz</</^/« 

(5) 

Sm{f/fm)~^ fm/ot<f<Otfm 



where = 2.7 x 10-'^^Hz-\ = 74 Hz, and a = 3.8 H], assumes particular detector 



parameters that accentuate the difference between Eqs. (j^) and (p. Since the ultimate 
shape of the noise curve is not yet well known, we feel that it is useful to calculate the 
attainable measurement accuracies for both of these shapes of the noise spectrum. We shall 
see below [cf. Tables II and III] that the flatter spectrum of the modified noise curve 
leads to a modest improvement in how accurately the binary's masses can be measured (for 
fixed signal-to-noise). 

B. Review of parameter estimation 

In this section we give a concise summary of those elements of signal processing that are 
necessary for parameter estimation. The basic concepts of detection and measurement have 
also been reviewed recently by Finn and by Krolak et al. in the specific context of 



laser interferometer gravitational wave measurements. In Appendix A we give give a more 
detailed treatment of parameter estimation, together with an extensive discussion of the 
ways in which the simplified linear formalism described in this section can break down: (i) 
when the signal-to-noise of the detection is low, and (ii) when our a priori knowledge of 
some of the binary parameters is not negligible compared to the information obtained from 
the measurement. 

We assume that an inspiraling binary gravitational wave has been observed; i.e, that 
the appropriate detection criterion has been met by the detector outputs. We now discuss 
how to determine the parameters of an inspiraling binary system that best fit the measured 
signal. The basic framework is illustrated in Fig. |I|. The set of all gravitational waveforms 
from two inspiraling bodies can be characterized by a relatively small number of parameters 
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(the distance to the source; the time of merger; five angles specifying the position of the 
source on the sky, the plane of the orbit, and the orbital phase at some given time; and 
the masses and spin angular momenta of the two bodies — fifteen parameters in all, assum- 
ing that the eccentricity of the orbit is negligible). We regard this set of waveforms as a 
fifteen- dimensional surface embedded in the vector space of all possible measured signals. 
In the absence of any noise, all measured signals from inspiraling binaries would lie on this 
submanifold; in practice, of course, the measured signal consisting of waveform plus noise is 
displaced off the submanifold. 

The statistical properties of the noise determine a natural inner product on the vector 
space of signals. Given two signals hi{t) and h2{t), we define {hi \ by p8| 



where hi and /12 are the Fourier transforms of hi and /i2- This definition is chosen so that 
the probability for the noise to have some realization no(t) is 

Hence if the actual incident waveform is h{t), then from Eq. ([I|) the probability of measuring 
a signal s in the detector output is proportional to g-C^'-'^l «-^)/2_ Correspondingly, given a 
measured signal s, the gravitational waveform h that "best fits" the data is the one that 
minimizes the quantity {s — h \ s — h)\ see Fig. |l]. 

It also follows from Eq. (P) that for any functions g{t) and k{t), the expectation value 
of {g\n){k\n) , for an ensemble of realizations of the detector noise is just {g\k). Hence 
the signal-to-noise (|]) of the detection will be approximately given by 

^\h]= ^^'^J , ={h\hYl\ (8) 
A^^ ^ rms {h\n) ^ ^ ' ^ ' 

The kernel w{t) of Wiener's optimal filter appearing in Eq. (0) is just the Fourier transform 

of l/Sn{f). 

For a given incident gravitational wave, different realizations of the noise will give rise 
to somewhat different best-fit parameters. However, for large S/N, the best-fit parameters 
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will have a Gaussian distribution centered on the correct values. Specifically, let ^* be the 
"correct" values of the parameters on which the waveforms depend, and let 9"^ + A^* be the 
best fit parameters in the presence of some realization of the noise. Then for large S/N , the 
parameter-estimation errors have the Gaussian probability distribution 

p{M''^)= Ne-"^^^^'-''''''. (9) 

Here is the so-called Fisher information matrix defined by 



F = 

^ I] — 



dh 



and Af = ydet(r/27r) is the appropriate normalization factor. It follows that the root- 
mean-square error in 6^ is 



^/{{AeTi = (11) 

where H = T~^. 

The above discussion applies to measurements made by a single detector. The (straight- 
forward) generalization to a network of detectors, which will be required in Sec. |rV|, is given 
in Appendix A. 

The above discussion also neglects the effects of any a priori constraints on the param- 
eters that may be available. The incorporation of such a priori information can have a 
significant effect on the predicted parameter-extraction accuracies (and also on the best-fit 
parameter values themselves). This is true not only for those parameters to which the con- 
straints apply, but also for the remaining parameters because of correlations. The effect 
is significant whenever, for some parameter, the a priori information is comparable with 
the information derived from the measured signal. Hence, a priori constraints are usually 
important whenever we include in an error-estimation analysis parameters which are weakly 
determined by the data. In Appendix A we derive a generalization of Eq. (^) [cf. Eq. ( |A45| ) 
below] which roughly incorporates the effect of a priori information. This generalization will 
be used in Sec. |111 lj| , where we consider the dependence of the inspiral waveform h on the 
spins of the two bodies. 
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C. The gravitational wave signal 
in the Newtonian approximation 



Inspiraling compact binaries can be described, to lowest order, as two Newtonian point 
particles whose orbital parameters evolve secularly due to gravitational radiation, where 
the gravitational waves and corresponding energy loss rate are given by the Newtonian 
quadrupole formula. That is, the orbital frequency Q at any instant is given by 

^ = 7^' (12) 

where M = Mi + M2 is the total mass of the system and r is the orbital separation. The 
inspiral rate, for circular orbits, is given by 

dr _ r dE _ 64 /iM^ 

where fi = M1M2/M is the reduced mass. Integrating Eq. (|13]) we obtain 

r=(^/iM2)'/'(t,-t)V^ (14) 

where tc is the "collision time" at which (formally) r ^ 0. Since the emitted gravitational 
waves are quadrupolar, their frequency / (cycles/sec) is equal to Q/n. The gravitational 
waves induce a measured strain h{t) at the detector which is given by (see, e.g., Ref. [|l^] ) 



Dr{t) 

where D is the distance to the source. The function Q and the angles 6, (p, ip, l (which 



describe the position and orientation of the binary) are defined in Sec. |V| below; they will 
not be needed in this section. In Eq. (|l^) we could have included the factor (384/5) ^^^tt^/^ 
in the definition of Q, but choose not to for later convenience. 

Because both the amplitude and frequency of the signal increase as t ^ tc, the signal is 
referred to as a "chirp." From Eqs. ( |T2|) and ([T3|) , the frequency evolves according to 

^ = ^^8/3 _^5/3 ^11/3^ (Ig) 
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where Ai = /i^/^M^/^ is the chirp mass parameter discussed in Sec. |. The phase of the 
waveform = 2tx f{t')dt' is 

0(t) = -2 [\M-\t, - 1)] + 0„ (17) 

where the constant of integration (p^ is defined by (/) ^ </)c as t — > te- 
la Eqs. (p!5D-(p!7D we have omitted the (obvious) time delay between signal emission and 
detection, and we have implicitly assumed that the detector and the binary's center-of-mass 
are at rest with respect to each other. The latter requires some explanation. If the detector 
and binary are in relative motion, the detected signal is Doppler-shifted with respect to the 
emitted signal. One cannot determine this Doppler-shift from the detected signal, since h{t), 
as defined by Eqs. ([T5| ) -([T6| ) , is invariant under the transformation 

{f,M,fi, r, D, t) ^ (//A, ^A, /iA, rA, D\, t\) . (18) 

Thus, strictly speaking, one can extract from the signal only the "Doppler-shifted" mass and 
distance parameters XA4, Xfi, and XD, where A is the Doppler-shift factor. This is not just a 
feature of our simplified, Newtonian waveform; it also holds for the true, general-relativistic 
waveforms, as can be seen on purely dimensional grounds and from the fact that general 
relativity does not define any preferred mass/length scales. 

Similarly, for binary sources at cosmological distances, the waves will depend on and 
reveal the redshifted masses 

M = {1 + z)Mtrne, = (1 + z)fitrne, (19) 

where z is the source's cosmological redshift, and also depend on and reveal its so-called 
luminosity distance Dl |]6|,|9|j26[| . Our measurement-accuracy analysis applies to these red- 
shifted masses and to the luminosity distance. The determination of the true masses for very 
distant binaries will require some method of estimating redshifts; see, e.g., Ref. 0. 
It is most convenient to work directly with the Fourier transform of h{t), 



oo 



oo 



h{f) = I e^^'f'h{t) dt, (20) 
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which is easily computed using the stationary phase approximation |T^. Given a function 
B(t) = A(t) cos (t>(t) , where d\nA/dt < d(f){t)/dt and d'^cp/dt'^ < {d(t)/dtf, the stationary 
phase approximation provides the following estimate of the Fourier transform B{f) for / > 0: 

B{f) ^ \A{t) (f )"'^'exp [^ {2nft - 0(/) - 7r/4)] . (21) 

In this equation, t is defined as the time at which d(f){t)/dt = 27cf, and (in a slight abuse of 
notation) 0(/) is defined as 0(t(/)). Using Eqs. ([TBI) and (1T7| ) we obtain 

t{f) = - 5(87r/)-«/3A<-5/3 

0(/) = 0, - 2 [87r^</]-^/^ (22) 

Hence from Eq. (|2l|) , the Fourier transform of the Newtonian waveform is 

h{f) = ^M'/'r'/'exp[^^{f)] (23) 

for / > 0, where the phase \I/(/) is 

vl>(/) = 2nft, -^c-\ + \{MMf)-"\ (24) 

Equation ( P^D for /i(/) is clearly invalid at very high frequencies, because the real inspiral 
will terminate at some finite orbital frequency. For BH-BH and BH-NS mergers, there will 
be a transition from inspiral to a final plunge near the location of the last stable circular 
orbit, which is roughly at r = 6M for non-spinning bodies The final plunge will last 



roughly one orbital period. [Neutron stars merging with rapidly spinning black holes may 
instead tidally disrupt, thereby shutting off the waves, outside the horizon [1^.] For NS-NS 
mergers, the two bodies will collide and coalesce at roughly r = 6M. Generally therefore 
the inspiral gravitational wave h{t) will "shut off" at roughly r = 6M, and correspondingly 
h{f) will shut off at roughly / = (G'^^^vrM)"^. We therefore "correct" the waveform (^3]) by 
setting h{f) = for / > {Q^/'^txMY^ . We note that when r > 6M, 



d(f)/dt 3{d(j)/dt)^ 55 'M' 
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so the stationary phase approximation should reproduce the Fourier transform of h{t) with 
good accuracy throughout the inspiral. Note that, as advertised in Sec. |, in the Newtonian 
approximation the signal ( p3D depends on Mi and M2 only through the chirp mass Ai. 

Using Eqs. and (p3D, we can see how the signal-to- noise squared accumulates as the 
frequency sweeps upwards: 



In Fig. ^ we plot the integrand d{S/NY/df = 4|/;,(/) p/S'„(/), using the advanced detector 
noise spectrum (^. The shape of this curve is universal once the noise spectrum is given: 
the masses, the relative angles, the distance to the source, etc. affect only the overall am- 
plitude. (This is strictly true only for the "Newtonian" signal, but will remain true to a 
good approximation when post-Newtonian effects are taken into account.) While 90% of the 
cycles come between 10 and 40 Hz, and while most of the energy is released in the last few 
orbits at / > 200 Hz, we find that ~ 60% of the total signal-to-noise squared accumulates 
between 40 and 100 Hz, the frequency band in which LIGO is most sensitive. 

We now evaluate the Fisher information matrix (10). For measurements using a single 



detector, there are only four parameters on which the Newtonian signal depends: an overall 
amplitude A = {Q/D)Ai^^^, and Ai, tc, and 0c- The derivatives of h(f) with respect to 
these parameters (for / > 0) are given by 

h, ^ = 2mfh, (27a) 



c 



d\nA ' dt 

^ = -zk #77 = -^{ST^Mf)-'/' h. (27b) 

From Eqs. (|27| ) and the noise spectrum (^, it is straightforward to evaluate the Fisher 
information matrix (^) and its inverse S*-' General expressions for the elements of Fjj 



using Newtonian waveforms, valid for any detector noise spectrum, are given in Ref . |M . We 



will not reproduce them here. However, for purposes of comparison to our post-Newtonian 



results in Sec. |ITT| , we list the rms errors A^, AAi, Ate and A0c for the case of low- mass 
(e.g., NS-NS) binaries, assuming the approximate waveform (|23|) and the detector noise 
spectrum (H): 

A(ln^) = 0.10(-^) (28a) 
Ate = 0.40 (^-^) msec, (28b) 
A0, = O.25(^-i^)rad (28c) 

A(1„A.) = 1.2 X 10- (^) (^) . (28d) 

For low-mass binaries, the fact that we "cut off" the waveform h{f) for / > (6^/^7rM)~^ has 
little effect on the rms errors (pHj), due to the sharp rise in Sn{f) at high frequency. The 
exact scaling of A(ln7Vl) as Ai^^^, and the fact that Ate and A0c are independent of M, 
strictly hold only when the cut-off is unimportant. For BH-BH binaries with S/N = 10, one 
has Ate = 0.60 msec, A0c = 0.32 rad, and A{\nM) = 1.3 x 10"^ (M/Mq)^^^. 

The rather phenomenal accuracy attainable for the chirp mass Ai is due to the large 
number N'cyc of cycles in the detectable portion of the gravitational waveform. We see from 
Eq. dT^ ) that Mcyc scales like Ai~^^^, so A(lnA1) is proportional l/A/'cyc, as one would 
expect. 

The rms errors ([28| ) apply to single-detector measurements. In practice, one will have 
a network of detectors, with different locations and orientations. For a network, A(lnA^), 
will be roughly given by Eq. (|28d|) , but with S/N replaced by the combined signal-to- noise 
p of the detector network, defined by Eq. (^ above. This is because independent estimates 
of Ai are obtained from each detector. The same argument does not apply to the rms errors 
in te, A, and (pc, because the gravitational waves will arrive at the different detectors at 
different times, and because detectors with different orientations measure different values of 
A and (j)e [cf. Sec. |rV| below]. 

We conclude this section by noting that from the measured value of the chirp mass Ai 
alone, one already obtains a lower limit on the larger of the individual masses, and upper 
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limits on the smaller mass and on the reduced mass. We adopt the convention that Mi > M2; 
i.e., Ml always refers to the larger of the two masses. Then it follows by definition that 

Ml > 1^1^ M, M2 < 2^1^ M, II < (29) 

However, if /i is unknown, then the mass ratio M1/M2 is unconstrained. The bounds ( pUj ) 
that follow from measuring M. may themselves be of astrophysical interest. For instance, if 
one determines using (^91) that Mi > 3Mq, then one may conclude that the heavier body is a 
black hole (assuming the redshift is small, cf. Eq. ([ToD above and associated discussion). Also, 
it has been suggested |§| that from LIGO/VIRGO measurements of NS-BH coalescences 
where the BH is rapidly spinning, it may be possible to constrain the neutron-star equation 
of state by measuring the frequency at which the NS's tidal disruption causes the waves to 
shut off. Knowledge of this tidal- disrupt ion frequency, coupled with an upper limit on the 
neutron star mass M2 determined from the inspiral waveform, would allow one to place an 
upper limit on the stiffness of the equation of state. 



III. POST-NEWTONIAN EFFECTS 
AND PARAMETER ESTIMATION 

We now extend the analysis of the previous section to include post-Newtonian effects. 
We continue to treat the bodies as point masses, since tidal interactions have a negligible 
effect. Also, for the moment we will neglect the effects of the bodies' spin angular momenta. 

The post-Newtonian approximation provides the most accurate description currently 
available of the gravitational radiation from inspiraling, stellar-mass binaries. Corrections of 
order M/r {P-^N corrections) to the lowest-order, Newtonian waveform (|15D were calculated 
almost twenty years ago by Wagoner and Will p2l. Calculations of the inspiral rate have 



recently been extended to P'^'^N order, for the case of non-spinning bodies, by Wiseman |33 



(after Cutler et al. |jl3[ and Poisson |Q had determined the form of the P^'^N correction for 
the case /x/M <^ 1). By ^^P^N order" we mean that corrections to the quadrupole-formula 
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radiation field and corresponding inspiral rate that are of order {M/rY have been taken into 
account, along with order {M/rY corrections to the non-radiative orbital equations which 
determine, e.g., the orbital frequency at a given separation. [There is no standard convention 
for "counting" post-Newtonian orders in calculations involving radiation; e.g., some authors 
refer to the lowest-order radiation field as P'^'^N. Our own terminology is motivated by 
the application considered here: since radiation reaction effects cause the inspiral, 0{M/r) 
corrections to the quadrupole formula accumulate secularly and have just as large an effect 
on the phase of the orbit as do 0{M/r) corrections to the orbital frequency at a given 
radius.] 

The post-Newtonian waveforms improve upon their Newtonian counterparts in three 
respects ||35|j36[| . First, they include contributions from higher-order multipoles of the stress- 



energy tensor (e.g., mass-octupole and current-quadrupole radiation in addition to the mass- 
quadrupole term), whose frequencies are different harmonics of the orbital frequency. Sec- 
ond, they include post-Newtonian corrections to the lowest-order expressions for the ampli- 
tude of each multipole component. And, most importantly for our purposes, post-Newtonian 
corrections to the energy E{r) and gravitational wave luminosity dE/dt{r) modify the in- 
spiral rate and thereby the accumulated orbital phase We can write h(t) schematically 
as 

h{t) = ^\j:hUt)e-''^A (30) 

L x,m ) 

where "3?" means "the real part of" , x indicates the term's post-Newtonian order, the integer 
m labels the different harmonics, and $(t) is the orbital phase. Each amplitude has the 
form 

hUt) ^ fmiMi/M,) QIXQ, ^, ^, (31) 

where r(t) is the orbital separation, is some function of the mass ratio, and is 
a function of the source's position on the sky and the orientation of the orbital plane. To 
connect with the notation of Sec. |T| and below, we note that the phase 0(t) of the quadrupole 
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part of the waveform is essentially twice the orbital phase: (f)(t) = 2^{t)+k, for some constant 
k that depends on the relative positions and orientations of the detector and the binary . 
Thus, the expansion for h{t) through P^'^N order is given by 

h{t) = 3fJ[ {hi + hl + hl-')e^''' + (/if + /i}-5)e** 

+ + /i^-5)e=^'* + hy"^ + /i^-^e^^*] (32) 

where $(t) has the post-Newtonian expansion 

$(t) = $0 ^ $1 + $1-5 + 0{M/rf. (33) 

In Eq. (^3]), $^ refers to the P^N order contribution to the orbital phase. As indicated by 
Eq. (^), the term ^^'^ vanishes identically, as do several omitted terms in Eq. (^21). The term 
f^o^2i<s>° jg j^g|- ^YiQ Newtonian, mass-quadrupole waveform given by Eq. (p!5[) , while the terms 
^o.5gi*o(t) ^^^^ ^o.5g3«$''(i) ^Yie lowest order current-quadrupole pieces of the waveform. 
The term h^'^ e^** is the so-called "hereditary" or "tail" term produced by the interaction 



of the outgoing wave with the binary's gravitational potential p^ , |38[l . The interested reader 
can find explicit expressions for the amplitudes through P^N order in Krolak [^. 

In Sec. I we argued that the waveform's accumulated phase $ contains most of the 
"information" that allows sensitive measurement of the masses of the bodies. Since this 
paper aims at only an approximate calculation of parameter-estimation accuracies, rather 
than use the full P^'^N waveform (^2|), we calculate the Fisher information matrix ([l0|) using 
the following "model" waveform: 

/i(t) = 5R|/i°e24*''+*^+*"l|. (34) 

That is, we include P^N and P^'^N corrections to the phase of the waveform, since these are 
decisive for extracting the mass and spin parameters of the binary, but we neglect the other 
post-Newtonian effects that are nominally of the same order. We expect that the values of 
AMi and AM2 calculated using Eq. ( |3^ will be a reasonable approximation to the error 
bars one would calculate using the true, general relativistic waveforms (assuming one had 
access to them). 
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There is another, practical, reason for the use of the truncated waveform (0). As 
explained in Sec. 0, we can simplify the error-estimation analysis by considering only single- 
detector measurements, and still obtain a reasonable estimate of the accuracies attainable 
for mass and spin measurements. However, as stated above, each of the amplitudes has 
a different dependence on the angles {9,ip,ip,L). These angles cannot be measured using 
one detector alone. The position of the source 6, (p is determined from differences in signal 
arrival times at (at least) three widely separated detectors P]. Moreover at least two of 
the detectors must have different orientations to obtain even a crude estimate of the angles 
ip, L (which describe the principal polarization axis of the wave and the angle between the 
line of sight and the normal to the orbital plane — see Sec. |I^ below). Thus, to make 
use of the extra information contained in the post-Newtonian terms that we are omitting 
in Eq. (0), a full detector network would have to be modeled. Hence, for simplicity, in 
our model waveform (|34D we omit all of the terms in Eq. (^) except for the largest one. 
[Although we do analyze a general network of detectors in Sec. fV| below, that analysis takes 
advantage of the fact that the phase-evolution information and the amplitude/polarization 
information in the measured waveforms are largely independent, and — complementary to 
this section's analysis — focuses on the amplitude/polarization information alone.] 



A. Parameter estimation neglecting spin effects 

In this section we estimate how well the masses Mi and M2 could be determined from 
the waveform, if we knew a priori (or a posteriori by some independent means), that both 
bodies had negligible spin. Note that this is different from the situation where the spins 
happen to be zero, but where we have no knowledge of this fact apart from the information 
contained in the gravitational waveform. 

In fact, it would not be justified to assume a priori that compact objects found in binaries 
have negligible spins. For one thing, the formation of close binaries generally involves a 
period of mass transfer, which would tend to spin up the accreting body. Observationally, 
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there are three known NS-NS binaries that will merge within a Hubble time; at the time 
of merger, the pulsars in these binaries will all be spinning at roughly 1 — 2% of their 
maximum possible angular velocities 0. (The spin rates at merger will be roughly a factor 
of two smaller than current values, due to magnetic dipole radiation.) We show in Sec. [Ill B 
below that allowing for spins of this magnitude increases the resulting error bars for mass 
measurements by roughly a factor of two, compared to the error bars obtained if spins are 
assumed to vanish. Nevertheless, we feel it is instructive to calculate the Fisher information 
matrix neglecting spin effects (i.e., assuming the spins are negligible a priori), both to 
illustrate the inclusion of post-Newtonian terms and to provide a basis for comparison with 
the results obtained when we include spins. 

We now briefly derive the P^'^N corrections to the phase of the waveform. Through 
pi-5jY order, the orbital frequency, energy, and energy-loss rates (for non-spinning bodies) 
are [021B3 
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where r is the orbital separation in DeDonder gauge (the standard gauge choice for post- 
Newtonian calculations), and t refers to time measured at infinity. 

Defining f = Q/ir, the frequency (in cycles/sec) of the quadrupolar part of the gravita- 
tional waves, we combine Eqs. (|35|)-(p7|) to obtain 



df/dt = 



1/3 



/743 llyu 
V336 ^ AM 



(nMf) 



2/3 



+47r (ttM/) + 0{TiMff'^ 



(38) 



In Eq. ( |5^ ) and below, we use (vrM/)^/^ as our post-Newtonian expansion parameter, instead 
of (M/r)^/^. We note that {tiM fY^^ equals [M/rY/"^ up to but not including terms of order 
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(M/r)^/^. This change of variables is advantageous because the frequency of the wave is a 
directly measurable, gauge- independent quantity (unlike the radius of the orbit). Equation 
(^) can be easily integrated to obtain t{f) and 0(/), where (p = it J fdt is the phase of the 
waveform. Defining x = (vrM/)^/^, we find that 

r 4 /743 

1 + \ X 

3 V336 4Mj 

(39) 
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where, as in Sec. ||, we define tc and 0c by t ^ tc and — > 0c as / — oo. 

Using Eqs. (|39D and (^) and the stationary phase approximation, we can repeat the 
analysis of Sec. ^ to obtain h{f). As before, we (crudely) model the end of the inspiral at 
r 6M by setting h(f) = for / > (G^/^ttM)"-*^. The stationary phase result then becomes 



j[f-7/6^i^ < / < (63/VM)-i 



(e^/VM)-! < /, 

where A = (Q/D) M^/^ and 



vI/(/) = 27r/tc - 0c - vr/4 + -(STrA^/)-^/^ 



(41) 
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Note that the post-Newtonian correction terms in square brackets in Eq. (^21) have their 
greatest effect on the phase of h{f) at low frequencies, because they are multiphed by the 
overall factor f~^^^. This may seem counterintuitive, since the post-Newtonian corrections 
to the inspiral rate are largest at small r, or high /; however the high-frequency portion of 
the waveform contains far fewer cycles, so the cumulative effect of PN corrections on the 
waveform's phase is smaller there. 

Our model waveform (^) for non-spinning bodies depends on five parameters: A, (pc M., 
fi, and tc- It is actually somewhat simpler to compute and interpret the Fisher information 
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matrix Tij in terms of the following modified parameters for which the rms errors are rescaled: 
In^, (pc, InA^, In/i, and /gtc, where /o is some fiducial frequency. With respect to these 
parameters, the derivatives of h{f) are |^ 

dh{f) 
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Using Eqs. ( ^Hj ) and the noise spectrum (^, we have numerically computed Tij, its inverse 
S*-', and the corresponding errors A(j)c = V^S^^^, etc. Since our model waveform includes 
post-Newtonian corrections to the phase but not to the amplitude, S*-' is block diagonal: 
Y^iriAj ^ Q fQj. j ^ InA^, In/i, or /otc- Hence, AA/A = (S'/iV)"^ while errors in A 
are uncorrelated with errors in the other parameters. Table | lists A0c, Ate, AM./ M., and 
Ayu//i for a range of values of Mi and M2. The results in Table | are for a single detector 
and are normalized to S/N = 10. For measurements by a detector network, the rms errors 
AJVijM., and Ayu/yU will be approximately those given Table |I|, but with S/N replaced by p, 
the combined signal-to-noise (|^) of the network. As explained in Sec. ||, this is because each 
detector provides almost-independent estimates of M. and /i. The result we particularly 
wish to draw attention to is: if spins can be treated as negligible then fi can typically be 
measured to ~ 1%, while Ai can be determined to ~ 0.01 — 0.1% 

Table |I| also lists the correlation coefficient = 5]'^^/(S'^'^ S^^')^/^, a dimension- 
less ratio indicating the degree to which errors in Ai and fi are correlated. The quantity 
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cmh is independent of S/N, and by definition satisfies CM^l ^ [^l?!]- We find that typ- 
ically \cm^i\ > 0.90, indicating that the errors in Ai and fi are strongly correlated. This 
strong correlation implies that there exists a linear combination of Ai and fi which can 
be determined much more accurately than either Ai or fi individually In particular, 



Ai^M- (S-^^/S^^) fij is smaller than AM by a factor of ~ (1 - c^^)"^/^. Indeed, the 
value of A (^A4 — (S-^'^/S^'^) jjj computed using our P^-^^ waveform ( Pl) is approximately 
the same as AA4 [cf. Eqs. (|2^) above] computed using the Newtonian waveform (^) pO| . 

How accurately can Mi and M2 be determined? While it is straightforward to answer this 
question when the mass ratio is large, we shall see that some care is required when Mi and 
M2 are comparable, since in this case the distribution of errors in Mi,M2 is non-Gaussian. 
Recall that we have adopted the convention that Mi > M2. Then we have 

1 
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Using Eq. (H), S^-^i^^ and S^^Ma ^^^^ expressed as linear combinations of S-^-^, S-^^, 
and S'*'^. However it is clear from Table | that in practice the S^'^ term will give the 
dominant contribution. Neglecting the terms proportional to S-^-^ and S-^^, we find that 
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_2/i(Mi - M2) 

For example, if Mi = IOM0 and M2 = 1.4M0, Eqs. (|5|) imply that AMi/Mi ^ 1.9A^/n 
and AM2/M2 ^ 1.4A/i//i. 

While the expressions ( ^5]) for S*"^^^^^ and J^^^^Ah should be adequate for estimating 
the distribution of errors when Mi ^ M2, these expressions unfortunately diverge when 
Ml = M2. This divergence is due to the fact that the Jacobian of the transformation 
(Mi,M2) {A4,fi) vanishes when Mi = M2. Of course, the rms mass measurement errors 



27 



do not actually become infinite. Rather, the linear approximation that one typically uses to 
estimate rms errors loses its validity. That is, the approximation that 

A r 9h ^ . dh ^ ^ dh ^ 

dh ^ ^ ^ dh ^ ^ ^ 

-AMi + — — AM2, (46) 



dMi dM2 

for variations Ah of a size determined by typical realizations of the noise, becomes inaccurate 
when Ml - M2 ^ 0, as dh/dMi + dh/dM2 ^ in this limit. 

To overcome this problem we proceed as explained in Sec. |A 7| below, and use the PDF for 
the best-fit values A^, /t of the parameters M., /i, which is a simply a Gaussian centered on 
the true parameters A^, jl. [Thus, we are considering so-called frequentist errors, cf. Sec. |A2 



below]. Let Mi and M2 be the corresponding best-fit values for the individual masses. 
Substituting into this PDF the transformation M = A4{Mi, M2) and fi = /i(Mi, M2) yields 
a non-Gaussian PDF for Mi,M2; from which we can calculate the 95% confidence limits 
for Ml and M2- The use of confidence limits is somewhat crude, in the sense that it leaves 
out much of the information contained in the PDF, but it is suitable for our purpose of 
determining roughly how accurately these quantities can be measured. Since AA^ is very 
small, for the purposes of this discussion we can assume Ai has been measured exactly. Let 
ft is the true value of the binary's reduced mass. Then with 95% confidence fi lies in the 
interval 

fl - 2A/i < fiK fl + 2A/i (47) 

where A/x = (S'*^)^/^ is determined from the variance-covariance matrix. Roughly speaking, 
a necessary condition for the distribution of Mi and M2 to be Gaussian is that /t + 2A/i < 
2~^/^ Ai (so that the yu's 95% confidence interval does not include the equal- mass case). 
From Eqs. (E4l) and (ff^) we obtain the following 95% confidence limits on Mi and M2: 



Mi(Al,/i + 2A/i) < Ml < Mi(Al,/i-2A/i) (48a) 
M2 ( Ai , /i - 2 A/i) < M2 < M2iM,jl + 2 A/i) (48b) 
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where the functions Mi and M2 are given by Eq. (|^) above. If /i + 2AyU is greater than the 
maximum allowed value of /i, then one should replace fl + 2AyU by 2~^^^Ai in Eqs. (^). For 
example, if M = 1.219M0, ft = O.TMq, and Afi/ft = 0.004 (the NS-NS case), then one can 
state with 95% confidence that Mi and M2 lie in the ranges IAMq < Mi < l.GSM© and 
1.2Mq < M2 < 1.4Mq. Thus Mi and M2 are determined with much less accuracy than n 
when the two masses are roughly equal. 

Figure |^ below illustrates the meaning of the confidence limits in Mi and M2 (in the 
context of more accurate calculations incorporating spin effects). 



B. Parameter estimation including spin effects 

We now present a rough calculation of the degree to which mass measurement accuracy 
is degraded when the spins cannot be assumed to be negligible. For the same reasons as 



in Sec. |III A| , we incorporate the effects of spins on the phase of the waveform, but neglect 
their effects on the waveform amplitude. 

Let 5*1 and S2 be the spin angular momenta of the two bodies, and let L be the total 
orbital angular momentum. We define the unit vector L by L = L/\L\. Then Kidder, 



Will, and Wiseman |^ have shown that, due to an "L ■ 5"' term in the two-body force 
law as well as spin corrections to the expressions for the system's mass-quadrupole and 
current-quadrupole moments, Eq. (p3) becomes modified at P-^'^N order as follows: 
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where again x = (nMfy^^, and where 
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Through P'^-^N order, the six components of Si and S2 affect the waveform's phase only via 
the particular combination ([50|) . (Of course, other combinations appear at higher order.) 
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We now discuss the magnitude of the correction due to p. For black holes, one has a 
strict upper limit on the magnitude of the spins: jiSjl < Mf. This is also roughly the upper 
limit for neutron stars, though the actual upper limit depends on the (uncertain) nuclear 
equation of state. We can therefore estimate the maximum size of (3 by considering the 
case where the spins are aligned with L, and where \Si\/Mf = |S'2|/M| = 1. In this case 
P = ■j^ — Y|(4/i/M). This maximum value Pmax is always within 10% of 8.5, regardless of 
the mass ratio. 

The P^'^N order equations of motion also contain "L x S" terms, which do not directly 
affect df /dt, but do so indirectly by causing the directions of L, Si and 5*2 to precess during 
the inspiral — essentially the Lense-Thirring effect. The equations describing the secular 



evolution of L, Si and 5*2 through P N order are [42 



dL 
'dt 



dSi 
~dt 
dS2 
dt 



^ aiSi + a2S2 

3 (^2 • L)Si + {Si ■ L)S2 



L 



X L 



aiLL + y2-l{S2-L)L 
a,LL + ^Si-^iSi-L)L 



X Si 
XS2, 



(51a) 
(51b) 
(51c) 



where ai = 2 + (3M2)/(2Mi), as = 2 + (3Mi)/(2M2), L = \L\ = ^i^fWr, and where, to this 
order, one can use the expression (|I4D for r{t). 

The precession of L, Si and §2 causes j3 to evolve; d[3/dt as calculated from Eqs. (0) does 
not vanish identically. Fortuitously, however, (3 is a/most conserved by Eqs. (0), in the fol- 
lowing sense. We integrated these equations numerically from / = lOHz to / = (G'^^^vrM)"^, 
for a wide variety of spin magnitudes, initial spin directions, and mass ratios; we found that 
(3 never deviates from its average value by more than ~ 0.25 (or ~ 0.03 /^^ax)- Moreover, 
the non-constant part of (3 is oscillatory, which further diminishes its integrated effect on 
the waveform's phase. These properties of the evolution of (3 are explored analytically and 
numerically in Appendix B. 

The near-constancy of (3 allows a considerable simplification of our model waveform: in 
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Eq. (^), we simply take /3 to be a constant. That is, we treat /3 as just another parameter, 
hke Ai and fi, on which the signal depends. The Fourier transform of our model waveform, 
including spin effects, is therefore given by 

A < / < (63/27rM)-i 

(63/27rM)-i < / 

where now 



(52) 



^if) = 2nft, 



(53) 



X 



20 /743 , ^ , 

1 + — ( — + TT^|x + (4/3- 167r) a;=^/2 



9 V336 4M. 

Now, for spinning bodies it is not really correct to treat the amplitude A = 
Q{9, 0, ip, l)D-^M^'^ as constant. The precession of the orbital plane described by Eqs. (|5TD 
causes the angles ip and l and to vary — and hence Q{6, ip, ip, l) to vary — throughout the in- 
spiral. Typically, the orbital plane precesses around the total angular momentum vector 
J = L + Si + S2 roughly 20 times during the observable portion of the inspiral. The re- 
sult is a sinusoidal modulation of the waveform envelope PJI3[], and the amplitude of the 
modulation can be large when l^il or 15*21 is comparable to \L\. Nevertheless, in the interest 
of simplifying the calculation, in our model waveform ( |52| ) we take A to be a constant. We 
discuss further below the implications of this simplification. 

The derivatives dh/dlnA, dh/d{fotc), and dh/dcpc of the signal ( ^21) are given by the 
same expressions as in Eqs. (^3]). The derivatives of h{f) with respect to InA^, In/x, and j3 
are: 



d\n M 

dhjf) 
91n/i 

dhjf) 
df3 



hi 



3i 
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{8nMf)-'/'h{f) 



55/x 



-3715 55/i 
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+ (247r-6/3)a;^/2 



3z {87rMf)-'/'{nMf)h{f). 



(54a) 

(54b) 
(54c) 
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Using Eqs. (0) we again compute the variance-covariance matrix S*-' for a range of values 
of Ml and M2. One can show that A^^, At^ AAi/Ai, and Afi/fi do not depend the value 
of p. The simple way to prove this is to make a change of variables from {A, (pc, tc,.M,fJ', P) 
to {A, (pc, tc,.M.,fJ', P') , where 

P' = (4/3 - levr) 7W^/V"^^^- (55) 

Since the waveform phase ^{f) [Eq. (|53|)] depends linearly on the Fisher information 

matrix calculated with respect to the new variables is independent of P' . This implies that 
the rms errors in the other parameters, and their correlation coefficients, are independent of 
the value of p. The values of A/?, cm/s and c^^ do depend on P, however. 

In Table |I|we list the rms errors A0c, At^, AM/Ai, Afi/ /x, and Ap for the same fiducial 
binaries that appear in Table |. In computing the results in Table H we use the model of 
the advanced detector noise curve given by Eq. (2.1). Since we are principally concerned 
with how our lack of knowledge of the bodies' spins affects how well we can determine the 
other parameters, we take the "true" value of P to be zero in all cases. As in Table |, the 
results in Table ^ are for a single detector and are normalized to S/N = 10; for a detector 
network, the rms errors AAi/Ai, and A^/fi and AP will be approximately those given in 
Table |Tl|, but with S/N replaced by the combined signal-to- noise ratio p. 

Summarizing the results of Table |I|, we find that AA^ / is roughly an order of mag- 
nitude larger than predicted by the Newtonian analysis of Sec. |I|, but still typically less 
than 0.1%. Thus, despite the "confusion" introduced by the extra parameters that enter at 
post-Newtonian order, we conclude that Ai can still be measured with remarkable accuracy. 
However, compared to the case where the bodies are assumed to have negligible spin a priori, 
we see that Afi has increased by a factor which ranges from 20 to 60! 

Table |T| also reveals the "reason" for this loss of accuracy: the correlation coefficient 
Cfj_j3 is extremely close to —1 Clearly the strong correlation is due to the fact that the 
frequency dependence of the "L ■ S" term in the expression ( p^D for the waveform phase 
\E'(/) is very similar to the frequency dependence of the other post-Newtonian terms in 
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Eq. (|53D. This strong correlation implies that there is a combination of /i and f3 which can 
be determined to much higher accuracy than /i itself [Q. Specifically, A (^/i — (S^^/S^^) (3^ 
is smaller than A/i by a factor of ~ (1 — c^^)~^/^, which is approximately 20 — 60 for the cases 
in Table |I|. Thus the combination fi — (S^^/S^^) f3 can be determined with approximately 
the same accuracy that one could achieve for /i, z/spin effects could be neglected [cf. Table 
H. Since both Ai and this particular combination of fi and (3 can be determined to high 
accuracy, the inspiral gravitational wave measurement essentially constrains the parameters 
to lie near a thin two-dimensional strip in {Ai,n,f3) space. This is illustrated in Fig. |^, for 
the case of a BH-NS binary. 

Up to this point, the formalism we have been using to calculate measurement accuracies 
neglects a priori constraints on the parameters, and thus implicitly assumes that f3 can take 
on arbitrary values. This assumption should be adequate as long as the 95% confidence 
intervals determined from S'^^ are well within the "allowed" range: \f3\ < j3max ~ 8.5. 
However we see from Table ^ that this criterion is not satisfied when both bodies are heavier 
than a few solar masses. For example, when Mi = M2 = IOMq we calculate Af3 = 19.5. 
We can (somewhat crudely) incorporate the restricted range of (3 into our formalism, as 
follows. We replace the a priori information \f3\ < /?max at hand by an assumed Gaussian 
distribution p^^\P) oc e~5(^/5)^ for /3. In Appendix A we derive an expression for the 
variance-covariance matrix which incorporates the effect of an (assumed Gaussian) a priori 
probability distribution for the signal parameters. We have used this result [Eq. ( |A45|) 
below] to re-evaluate the variance-covariance matrix for the two high-mass binaries shown 
in Table |I[ (Taking the restricted range of f3 into account makes little difference to the 
other cases in Table y.) These re-evaluated results are marked in Table ^ with a dagger 
(t). Again, the rms errors listed are for S/N = 10; note however that since p^^\P) is fixed, 
the rms errors no longer scale simply as {S/N)~^. We see that taking the restricted range 
of (3 into account leads to the improved estimate Afi/fi^ 50% in both the high- mass cases. 

We mentioned above that in the three known short-period NS-NS binaries, the radio 
pulsars will, at the time of merger, all have spin angular momenta that are ^ 2% of their 
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maximum possible values. We feel it is an interesting exercise to calculate what measurement 
accuracies could be attained if we knew that NS's in nature were slowly spinning in general, 
e.g., if we knew a priori that (3 < .02j3max for NS-NS mergers. Repeating the procedure 
used above, we take p^^\P) oc e~^(/^/o-i)^^ and we use Eq. ( [A45| ) to calculate the variance- 
covariance matrix for the NS-NS case, for S/N = 10. We find AfJ^/fJ^ ~ 0.9%, which is 
roughly twice the value obtained in Sec. |III A| , where spin effects were taken to be completely 



negligible. 

We turn again to the question of how accurately the individual masses can be measured. 
The procedure for calculating AMi and AM2 in terms of AA^ and Afi is of course the same 



as described in Sec. pi A] . Thus for the BH-NS case, using the fact that Ayu/yU ~ 15%, we 
find from Eqs. (H) that AMi/Mi ^ 30% and AM2/M2 ^ 20%. In Sec. we explained 
that the distribution of errors in Mi and M2 will be non- Gaussian if /i + 2 A fx > 2'^/^ M. 
By this criterion, if A/i//i ^ 15%, then we can reliably estimate AMi and AM2 by using 
Eqs. (H) only if M1/M2 > 5.5 

Again, even when the Gaussian approximation is invalid, one can still use Eqs. (^81) to 
place 95% confidence limits on Mi and M2. Consider again the NS-NS case, which we looked 



at in this context in Sec. |111 A|, with the true values of the masses being Mi = M2 = 1.4M0. 



Then Ai = 1.22Mq, and, using the 2(7 error bar indicated by Table |l| (for S/N = 10) 
we see that, 95% of the time, the observers would measure fi to be between O.56M0 and 
O.7OM0. Correspondingly, the measured values of Mi and M2 would lie in the ranges 
I.4M0 < Ml < 3.2M0 and O.7M0 < M2 < 1.4M0. Thus, in the NS-NS case, measuring 
fi to within 20% means determining the individual masses only to within a factor of ~ 2. 
The constraints obtained on Mi and M2, for this case and the BH-NS case, are illustrated 
in Fig. |. 

Finally, we repeat these calculations using the flatter spectrum instead of (||) as our 
model of the advanced detector noise. The results are shown in table We see that the 
main conclusions which we drew from table ^ are unchanged, but that (for fixed signal- 
to-noise) the relative errors AA4/A4 and Afi/fi are a factor of ~ 1.5 times smaller with 
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the flatter noise spectrum (|). This is presumably due to the fact that noise spectrum 
exhibits better sensitivity at low frequencies, where most of the gravitational wave cycles 
(and hence most of the sensitivity) come from. 



C. Caveats and future work 

Since the results in Tables |-|T| were obtained using several approximations and sim- 
plifying assumptions, we feel that it is useful to collect the most important of these in one 
place. They are as follows: 

First, we restricted attention to statistical errors arising from detector noise. In practice, 
theoretical template waveforms will be quite difficult to compute accurately ||8|jr3|l. Hence 
some systematic error may also arise from fitting the data to imperfect template waveforms. 
Currently a large effort is underway in the relativity community to calculate templates suffi- 
ciently accurately that these systematic errors will be at most comparable to the statistical 
errors that we have obtained — at least P'^N and possibly higher order templates will be 
required. We note that template inaccuracies, while giving rise to important systematic 
errors in parameter extraction, will not significantly diminish our ability to detect the waves 

i- 

Second, we assumed the "advanced LIGO" noise-curve shape, for which we have used 
two estimates: Eq. @ and Eq. (|^). As emphasized above, these are only rough estimates 
of the spectral shape that the LIGO/VIRGO detectors will actually achieve. We have seen, 
however, that our results do not depend very sensitively on the exact shape of the noise 
spectrum. 

Third, we have used the approximate, linearized error-estimation formalism described in 
Sec. |11B| and Appendix A; the rms errors so calculated are guaranteed to be accurate only in 



the limit that the errors are small. When the errors are so large that the linearized approach 
is invalid, then our approach will probably generally underestimate the true variances. To 
avoid the limitations of the linearized error analysis, we are currently performing a Monte- 
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Carlo simulation of the parameter extraction process. 

Fourth, we calculated the variance-covariance matrix S*-' using the simplified model 
waveform (^2]), which is qualitatively inaccurate in a number of respects. In particular, our 
model waveform depends on the spins of the two bodies only through a single parameter, 
(3. We have neglected the spin-induced precession of the orbital plane, which also arises at 
pi-SjY order, and we have neglected the effect of the spin-spin coupling on inspiral rate, which 
arises at P^N order. We have also neglected higher-order multipole radiation (except insofar 
as the energy carried away by the higher multipoles affects the inspiral rate), and have only 
crudely modeled the cut-off of the waveform during the bodies' final tidal-disruption, plunge, 
or coalescence. [In particular, we have made no attempt to model the spin dependence of 
the cut-off.] 

It is unclear to us whether the inadequacies of our model waveform have led us to 
underestimate or overestimate parameter-extraction accuracies. On the one hand, the P^N 
spin-spin interaction term that we have neglected would, if included, inevitably lead to 
some degradation of parameter-extraction accuracy (as always happens when there are more 
parameters to fit for). On the other hand, it seems clear that the effect on A/i of adding the 
spin-spin term to the waveform will be far less dramatic than the inclusion of the spin-orbit 
term /?, for two reasons: (i) as shown by Kidder et al. for the BH-NS and NS-NS cases. 



the effect of the spin-spin term on the accumulated phase of the waveform is a factor of at 
least 20 smaller than the effect of the spin-orbit term, and (ii) the correlation coefficient \c^f^\ 
is so close to 1 because the frequency dependences of the /x and {3 terms in the waveform 
phase '^{f) are so similar; the /i and spin-spin terms are less similar in their frequency 
dependence. 

Finally, by neglecting higher-order multipoles, spin-precession effects, and the details 
of the final plunge, we have effectively thrown away information that would be contained 
in the true waveform. In a more complete analysis, this "additional" information could 
perhaps decrease measurement uncertainties. In particular, if in some cases the spin-related 
modulation of the waveform carries substantial information about /?, then it is clear from 
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Fig. ^ that A/i and consequently AMi and AM2 could be reduced by large factors. This is 
an important possibility which we are currently investigating. 



IV. ACCURACY OF DISTANCE MEASUREMENTS 



A. Overview 



In the previous sections we have investigated how accurately the masses of the inspiraling 
compact objects can be measured from the phase evolution of the detected gravitational 
waveforms. The other interesting parameters that are measurable from the outputs of a 
network of detectors are the distance D to the source, and its position on the sky. These 
parameters will be encoded in the amplitudes, phases and arrival times of the signals ha{t) 
read out from the detectors. At least three geographically separated detectors will be needed 
in order to determine the distance |0. We start by describing, in detail, the dependence of 
the signals ha{t) on the binary's distance and sky location. 

Let Xq be the position and be the polarization tensor of the ath detector in a detector 
network. By polarization tensor we mean that tensor for which the detector's output is 
given in terms of the waves' transverse traceless strain tensor h(x, t) by 



Here the colon denotes a double contraction. If the arms of the detector are in the directions 
of the unit vectors / and m, then d„ = ( / Cg) / — m ® m)/2 [^. 

We introduce a spherical polar coordinate system (6*, ip) centered at the Earth so that 
the axis = is the Earth's axis of rotation. The angle ip is longitude and 71/2 — 9 is North 
latitude for 9 < 7i/2. Let n = (sin 9 cos ip, sin 9 sin ip, cos 9) be the unit vector in the direction 
{9, p), and let e+ and e^^ be a basis for the transverse traceless tensors perpendicular to 
n. If we demand in the usual way that : = 25^^ , for A,B = +,x, then this basis is 
unique up to rotations of the form 



hait) = da : h(xa,t). 



(56) 




(57) 
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The quantities 

F,^(n) = e;^ : d,, (58) 
for A = +,x, are the so-called detector beam-pattern functions for the ath detector [O. 



Consider a coalescing binary source in the direction n. As in Sec. [Ill .Bl, let L denote the 



unit vector in the direction of the binary's orbital angular momentum, and let 

V = cos 6 = L ■ n, (59) 

so that 6 is the inclination angle of the orbit to the line of sight. As seen from the Earth, 
the orbit looks elliptical, and the principal axes of the ellipse give a preferred polarization 
basis e'+, e'^ for the waves. Specifically, we define 

, n X L , ^ 

< = 7. ^, (60) 

||n X L|| 

< = JT^r (61) 

^ ||n X e'^,|| 

where the minus sign is inserted to accord with standard conventions; the waves propagate 
in the direction — n. The preferred basis is e'+ = (g) — e^, e'^ = e'^ ®e'y + e'y^ e'^. 
In terms of this basis, the waves' strain tensor is h{t) = h^{t)e'^ + hy.{t)e'^ , where in 
the quadrupole-moment approximation the waveforms h^{t) and hy.{t) are as given in, e.g., 
Ref. [|12|. Taking the Fourier transform we find 



hA{f) = XA{v)hoif), (62) 
where x+(f) = (l + f^)/2, Xx{v) = —iv, and 



for / > 0. The phase ^{f) is the same as previously given in Eq. (|5^), and depends only 
on the parameters A^, /i, /5, tc, and 0c- 

If we fix a polarization basis e"*", e^, then we have 
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JA 



(64) 



for some polarization angle ip, where R\ is the rotation matrix 



^ cos(2^) sin(2^) ^ 



V 



(65) 



— sm{2ip) cos{2ijj) 

The conventional definitions of e+, and the corresponding definition of ip for a single 



detector are given in Refs. |ll2| , ^ . A network of several detectors, however, determines a 
different preferred basis e"^, (see below), so for the moment we allow the basis to be 
arbitrary and define ip via Eq. (|64D. By combining Eqs. (p6D, (0), (0), and ( |64D we obtain 
the signal read out from the ath detector: 

haif) = R%{2^P) XA{v)Ff{n) e^-^-^ ho{f), (66) 

where = — n ■ The first three factors in Eq. ( |66D taken together are proportional 

to the quantity Q{9,ip,'ip, t) that appears in Eqs. (|I5]) and ([231). 

Now the overall amplitude A of the signal at one detector can be measured to an accuracy 
(cf. Sec. 0) 

where pa = {ha \ ha)^^'^ is the signal to noise ratio (SNR) measured at that detector. Since 



A oc 



Q{9,ip,ip,L) 
D 



(68) 



we expect the accuracy of distance measurements to be very roughly AD / D 1/p, where 



X] Pa is the SNR ( |A27|) , giving an accuracy of ~ 10% for typical detected signals. 



However from Eq. (^) the signal amplitudes are also strongly affected by the angles 9, if, 
ip, and most importantly the inclination angle l. Hence, there will be correlations between 
the measured values of D and of these angles, and the accuracy of distance measurement 
will be reduced relative to the above naive estimate based on Eq. (|67D. 

It is straightforward in principle to calculate the effect of all the correlations by calcu- 
lating the Fisher information matrix (|lOD from the waveform (^6]) for all of the variables 
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D, n, and together with the variables A^, 0c, ^c, and /5 discussed in Sees. |T| and 
pn] . An analysis of this sort, but without including the post-Newtonian parameters /i and 
/?, has been carried out by Jaranowski and Krolak ||I5[, who numerically calculate the rms 
error AD for various different values of the angular variables. They use the three detector 
network consisting of the two LIGO detectors and the VIRGO detector in Pisa, Italy, with 
their planned orientations. Although these authors do not take into account post-Newtonian 
effects, it seems likely, for reasons which we discuss below in Sec. [IV Q and Appendix C, that 
their results for /\D will not be sensitive to this restriction. Similar numerical calculations 
have been carried out by Markovic who assumed the same network of detectors. He 
identified a useful approximation for calculating AD/D^ based on identifying those vari- 
ables with which the distance measurement is most strongly correlated, and neglecting the 
effect of the much smaller correlations with the other variables. 

In Sec. [IV Q below we present an analytic calculation of AD/D which simplifies the 
treatments given in Refs. P,p!5|. Because the rms error AD depends on several angular 
variables, it is difficult to explore its behavior over the whole parameter space using numerical 
calculations of the type in Refs. P,p3|. Here, by using Markovic's approximation, we derive 
an approximate analytic expression for AD, which is valid for any network of detectors. 

We also extend the analysis of Refs. P,|l5i in the following two respects. First, we 
parametrize the dependence of the result on the positions and orientations of all of the 
detectors in the following useful way. We show that, for a given position 9, Lp on the sky, 
the detector network parameters influence AD only through (i) the selection of a preferred 
polarization basis (e+,e^) [or equivalently a preferred polarization angle V'(n), cf. Eq. ( [57[) 
above], and (ii) two quantities aoip) and £z)(n), where we call aoij^) the amplitude sensi- 
tivity and 1 — £^(11) the polarization sensitivity [^. We discuss these "network sensitivity 
functions" in detail in Sec. [IV Ij| below. They are defined in such a way that the total 
signal-to-noise ratio squared ( [A27[ ) of a detected signal coming from the direction n with 
polarization ip is of the form [cf. Eq. (RTI) below] 
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oc crD(n) [1 + £z)(n) cos(4^ + const) f{v) ] , (69) 
where the function 

(1 -f2)2 

is independent of n and ip. The values of ct/j and e^) are, as an example, ctd = 1 and = 
for the case of two detectors at the same location, rotated with respect to each other by 45°, 
and for vertically incident waves. In Figs. ^ and ^ below we show plots of these quantities 
as functions of the angles 6 and (p, for the 3-detector, LIGO/VIRGO network. 

Second, we extend in Sec. [IV D| the analysis beyond the linear, Gaussian approximation 



outlined in Appendix A, which is normally used to estimate the rms errors. We do this 
by calculating the exact (within the Markovic approximation), non-Gaussian probability 
distribution for the distance D which incorporates both our a priori knowledge and the 
information obtained from a gravitational wave measurement. This extension becomes im- 
portant in two different regimes. The first regime is when v = cost 1, corresponding 
to binaries that we perceive to be almost face-on. In the limit f — 1, the value of AD 
predicted by the linear approximation becomes infinite. As shown by Markovic this is 
because two of the signal parameters become degenerate (i.e., the derivatives dh/d6^ become 
linearly dependent) as f ^ 1. Markovic gave rough estimates of the effect of this breakdown 
of the linear formalism on the predicted value of AD; the effect is not treated in the exact 
numeric calculations of Ref. |jl5[. Here, using the non-Gaussian distribution for D, we ob- 



tain an improved approximation to AD near the points of degeneracy. The second regime 
where our non-Gaussian extension of the error-estimation method is important is the limit 
of low signal-to-noise, and correspondingly of large relative errors in the measured binary 
parameters. Since the Fisher matrix method of calculating the rms errors in the measured 
parameters gives essentially the leading order term in an expansion in powers of {S/N)~^, 
this method will be inaccurate at low values of S/N. By using an approximation which takes 
into account the dominant effects that are non-linear in (S/N)^^, we numerically estimate 
AD for different values of the parameters n, v, and ip- We show that the linear estimates 
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for AD are typically off by factors ^ 2, even for signal-to-noise ratios of more than twice the 
threshold value for detection, due in some cases to large non-Gaussian tails in the PDF for 
D. Thus, effects that are non-linear in {S/N)~^ are often not a small correction for typical 
detected signals. 

Finally, in Sec. |1V ii| we apply our non-linear error estimation method to calculate the 
distribution of measurement accuracies for the LIGO/VIRGO network, using a Monte-Carlo 
simulation. We estimate that ~ 8 % of the distance measurements will be accurate to < 15 %, 
and ~ 60 % to < 30 % (see Fig. [l| below). 

Our analyses are applicable to binaries at cosmological distances, provided we interpret 
D as the luminosity distance to the source, and M. as + times the true chirp mass, where 
z is the source's redshift [|,mi4|. However, a potentially important effect that we neglect is 
the spin-induced modulation of the signal amplitudes discussed in Sec. [Ill B| and Ref. 



Hence, our results for /S.D should be regarded as rough estimates (and probably also as lower 
limits, since it seems most likely that including spin effects in the computation will increase 
AD). However, the tools we develop below will be useful in future, more complete analyses 
of distance measurement accuracies. 

We use throughout this section the notations of Appendix A. 

B. The network functions ^"^(n) and eD(n) 



The overall SNR (|A27|) and the Fisher information matrix (|A31|) are determined by 



inner products involving the signal h(t) and its derivatives dh/d6^ with respect to the signal 
parameters 6\ We now show that a large class of these inner products depends on the 
network properties [i.e., the detector positions Xq and polarization tensors da] only through 
the two functions of sky location, a£){n) and ^^(n). We start by defining the complex 
amplitudes 

AB^RV2^)XAiv)e-'^^/D (71) 
which are intrinsic to the incident waves, and the detector amplitudes 
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A= E •^BF^i.n) (72) 

-B=+,x 

which characterize the signals seen at the various detectors. In terms of these quantities, 



the signal ( |66D can be written as 

ha{f)=Aae'-'^^f~kU). (73) 

where k{f) = De^'^^ho^f) is independent of D and 0c- The inner product of two signals h 
and h' written in this way, with amplitude parameters Aa and A!j^, is given by Eqs. (|63[), 
(0) and 

(h I h') = 3? [^M'fe k"^] {k\k). (74) 

Here the positive definite Hermitian matrix k"-^ is 
.6 ^ S^dfri^ [SMY'f e^-/(----) 

lo^dff-yysM) ' 

and Sn{f) in the denominator is the average of the spectral noise densities in all the detectors. 
If the detectors are all identical, and correlated sources of noise [represented by the off- 
diagonal elements of S„(/)] are unimportant, then k*^* is just 5"*. 

In terms of the wave amplitudes Aa, the inner product (ffiD is, from Eq. ([7^), 

(h I h') = 3? [A*aA'b e^^] (^ I ^) , (76) 

where the matrix is given by 

©^^(n) = E^a^(n)n^(n)K"^ (77) 

a,b 

We see that all inner products of the type (^) depend on the network parameters only 
through the 2x2 Hermitian matrix 0. Two key simplifications now arise. First, correlated 
sources of noise will presumably be limited to pairs of detectors at the same detector site, 
so that the detector-network noise matrix ( |X3| ) will have a block-diagonal form with each 
block corresponding to a detector site. If the detectors at each site are all oriented the same 
way, as is likely, then the product of beam pattern functions appearing in Eq. (|77|) will be 

43 



constant over each block in the indices a, b that corresponds to a non-zero subblock of the 
matrix S„(/). Hence, from Eqs. ( fTSD and (|77D, we see that the imaginary part of will 
vanish. Second, if we change the basis e+, by a transformation of the form (pTj), which 
amounts to redefining the polarization angle by 



^ = + Alp, 



(78) 



then will transform according to R(2A'?/^) ■ ■ R(— 2A'^). For fixed n, we can by 
choosing Aip suitably make diagonal, and so be of the form 



6 = an 



+ eD 



(79) 



l-SD 

where < e/j < 1. This defines the network functions cr/)(n) and eoi'n). The required value 
of Atp = Atpln) is given by 

2e+x 



tan(4A^) 



(80) 



0++ + exx 

The combined SNR ( [A27] ) can be determined in terms of these network functions by 
combining Eqs. ([FTD, (|7|) and ( |79D to give 



co(f) + enin) ci{v) cos(4^) 



(81) 



where Co{v) = (1 + v'^^/A + f ^, Ci{v) = (1 + v'^^/A — v'^ = {1 — f^)^/4, and ip is given by 
Eqs. (|7|) and (H). The quantity 

\ho{f)\' 



pi = D-\k\k) = (/io|/io)=4 



SnU) 



df 



(82) 



appearing in Eq. (|8T|) is the SNR that would apply to one detector if a face-on {v = 1) 
binary were directly overhead. From Eqs. (^) and (0), we find po = i^o/D, where for the 
noise spectrum (H) the distance vq is 



ro = 6.5 Gpc 



5/6 







-2/3 



70 Hz 



So 



3 X 10-48 sec 



-1/2 



(83) 
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The fiducial values of the detector parameters 5*0 and /o used here are those appropriate for 
the advanced LIGO detectors 0, cf. Sec. y above. Note that the dependence of the SNR 
(pT|) on the polarization angle ip vanishes when the binary is perceived to be face-on {v = 1), 
as we would expect physically due to rotational invariance about the line of sight. In the 
opposite limit of edge-on binaries {v —* 0), the incident waves are highly linearly polarized, 
and the SNR typically depends strongly on ip, varying by factors of ~ 10 or more as ip is 
varied [see Fig. ^ below]. 

We now derive simple formulae for the functions ao and Ed- From Eq. (|79D , it is clear 
that a£) is just half of the trace of the matrix 0, which is invariant under rotations. Using 
Eqs. (|5^) and (|77D gives 

= |e (e^)., (e;^)..} |e (d.)., (dj^j • (84) 
If we denote the first term in curly brackets by 5'^^^^, then it is straightforward to show that 
Sijki = -SijSki + {SikSji + 5u5jk) + {SijUkTii + SkiniUj) 

+ninjnkni. (85) 

This yields for the amplitude sensitivity function the formula 

anin) = i/t"^ [2d„ : d, - 4n ■ (d„ ■ d,) ■ n 

+ (n-d, ■n)(n-d5-n)] (86) 

where we have used the property Tr da = 0. 

It is similarly straightforward to evaluate the polarization sensitivity 1 — eoln). We 
introduce the notation 

(da I db)^ = {da)ij Sijki {db)ki (87) 
= 2d^d^-(Trd;[)(Trd^), (88) 
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where denotes the projection {Sik — nink){Sji — rijUi) {da)ki of d^ perpendicular to n. 
Then, using the relation from Eq. (0) that Tr 0^ = 2cr|,(l + ef,), and Eqs. (|8D and 
gives 

eoirvf = y\t-, E (da I d,)„ (d, I d,)„ - 1. (89) 

^ V / abed 

We now evaluate an and eo for the LIGO/VIRGO detector network. Let e^; = n, 
and be the usual basis of orthonormal vectors. Then for a detector at position 6, ip on 
the Earth's surface, such that the angle measured anticlockwise from the local eastwards 
directed meridian to the bisector of the detector arms is a, the polarization tensor is 

d = - sin(2Q;)(eg ® - ® e^)/2 

+ cos(2a) (e^- ® + ® eg)/2. (90) 

The values of {9, if, a) for the various detectors are (59.4°, —90.8°, 243°) for the LIGO de- 
tector in Hanford, Washington, (43.5°, —119.4°, 171°) for the LIGO detector in Livingston, 
Louisiana, and (46.4°, 10.25°, 117°) for the VIRGO detector in Pisa, Italy |T§. We assume 



that the detectors at all three sites are identical and that noise sources are uncorrelated, so 
that from Eq. ([7^), = 5"^. The resulting plots of and I — Sd are shown in Figs. ^ and 
i 



C. The Markovic approximation 

We now explain the approximation method used by Markovic 0, which we modify 
slightly below. We start by considering the accuracy An with which a given source can be lo- 
cated on the sky. The location n will be largely determined by "time of flight" measurements 
between the various detectors, i.e., measurements of the quantities — = — n ■ (x^ — X5) 
in Eq. (pBf ) @JT^. Hence, the variables n and will be strongly correlated, and An will be 
largely determined by the ratio of the timing accuracies A(ra — r;,) to the light travel times 
between the various detectors. Schutz ||^ has estimated the resulting angular resolution to 
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be ~ 1 square degree for typical detected signals, which is roughly in agreement with the 



recent detailed coalescing binary calculations of Jaranowski and Krolak |T^. It is also in 
rough agreement with numerical simulations of Giirsel and Tinto , which were carried out 
in the context of arbitrary bursts of gravitational waves. Hence, we see from Eq. (p6| ) that 
typical variations in n will give rise to variations in the measured value of D that are small 
compared to AD. Thus, the correlations between D and n should be small, and to a good 
approximation we can treat n as known when calculating /S.D |^. 

In the approximation that n is constant, we can divide the remaining parameters into 
two groups. The first consists of the four "amplitude" parameters D, f , if)^ and 0c, which 
determine the two complex amplitudes A+ and via Eq. (^TJ). The second group of pa- 
rameters consists of M., tc, /i, together with some spin parameters, which enter only in the 
phase of the Fourier transform of the signal, and control the evolution in time of the 
phase of the waveform In Appendix C we show that the second group of parameters 



decouples from the first to linear order in 1/p, and in the constant n approximation. More 
precisely: if one calculates the Fisher matrix ( [A31D for all of the variables except n, inverts 



it to obtain the covariance matrix E*-', and takes the 4x4 subblock of S*-' corresponding 
to the amplitude group of parameters, then the result is the same as if one computes the 
Fisher matrix for just the four amplitude parameters alone, and then inverts that. Heuris- 
tically what this means is that the effect of the correlations between [D, v, tp) and all of the 
parameters (j)c,tc,A4, etc., can be computed by considering the correlations with just one 



phase variable, namely 0c; the orbital phase at coalescence [49 



We now calculate the Fisher information matrix (|A3 1| ) for the four amplitude param- 



eters D, V, ip and 0c, and for an arbitrary detector network, as this should yield a good 
approximation to AD. The approximation that was used in Ref. |^ was in fact to consider 
only V, and tp; below we find [cf. Eqs. ( p2| ) and (p^)] that the fractional corrections due to 



also including 0c are of order eD(n) sin(4^/'), where ijj is given by Eqs. (|78D and (^). Since 
< Sd < 1 always, the fractional corrections are always ^ 1 [Q. 

From Eq. (|A31|) , it is clear that the Fisher matrices calculated using two different sets 
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of variables are simply related by transforming with the Jacobian matrix of the variable 
transformation. Hence, we can use any convenient set of variables to evaluate Tij and S*-^, 
and afterwards transform to the physical variables of interest. We define the variables a and 

( \ 



1 



(91) 



The waveform (|66D depends linearly on these variables, which simplifies the computation. 
Using Eqs. - (|7|), (|7|), (|7|), (|8|), (|A3TD and the relation po = r^l^, we obtain 



rl (T^(n) [Fij- +eB(n) G 



u J ' 



(92) 



where the variables are = (■?/', a, P, (j)c)- Defining cq = + Ci = — a^, C4 = cos(4-?/'), 
and S4 = sin(4'?/;), the matrices F and G are given by 



4co 

1 



-4a/3 





1 




-Aa/3 


Co 



(93) 



and 



G 



(94) 



-4ciC4 2as4 —2[3si 
2a;S4 — C4 —Ps^ 
— 2/3s4 C4 asi 
— /3s4 as4 C1C4 
Inverting the matrix ( ^2]) and taking the 2x2 subblock corresponding to the variables a, /3, 
we find with the help of Mathematica that 



1 



1 + £dC4 

1 — EoCi 



(95) 



Finally, transforming this with the Jacobian of the transformation ( PT| ) and taking the (D, D) 
element of the resulting matrix yields 
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Q T~)4 

AD2 = ^T(n,t;,^)2, (96) 
rid ro 



where 



(1 + v^) - Ed cos(4V^) (1 - v^) 



and rirf is the number of detectors. The dimensionless function T satisfies 

T(n,t;,^)>l, (98) 



since from Eqs. (p6D, (|88|), and ( [79D it follows that (Td < and < Ed < I always, for 

any detector network. 

Equation (p6| ) is the main result of this subsection. We now discuss its properties and 
range of applicability. It clearly breaks down and overestimates AD when v —>■ 1. As 
shown by Markovic, this is because dh/dD oc dh/dv at i; = 1, so that the linear error- 
estimation method breaks down. However, it will underestimate the true measurement error 
for sufficiently small values of the SNR p, because of the inadequacy of the linear error- 
estimation formalism in this regime (cf. Appendix A). In Sec. [IV D| below we numerically 
calculate more accurate values of AD/D, and show that even for small values of v, and 
even for relatively large values of p (e.g., p ^ 20, more than twice the threshold), the results 
predicted by the formula (|96D can be off by factors ^ 2. 

Hence, the formula (|9^) is of only limited applicability. Its main virtue is that it allows 
one to understand qualitatively how the distance measurement accuracy is influenced by the 
parameters ctd, £d, 'ip, and (to a more limited extent) v; and thereby by using Figs. |^ and 
^ how it varies with sky location n. We now discuss the dependence of AD /D on these 
parameters. 

As the polarization angle ip is varied, it can be seen that 

Tj„in(n, v) < T(n, f , ip) < T^^^{n, v), (99) 

where Tmin and Tmax are given by substituting cos(4'?/;) = ±1 in Eq. (P7|). As an illustration, 
Figs. and || show Tmin and Tmax as functions of n at f = 1/ -\/2- It can be seen that the 
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distance measurement accuracy can vary over the sky by factors of order ~ 20, for binaries 
at a fixed distance and with fixed inchnation angle. The reason for this strong variation of 
more than an order of magnitude is easy to understand. A key feature of the resuh ( pTj) 
is the factor of 1/(1 — ej^), which diverges in the hmit ^ 1. This divergence is not an 
artifact of our approximate, hnear error-estimation method (unhke the divergence in T at 
f — ^ 1). The physical reason for the divergence as £:d — » 1 is that for directions n such that 
1 — ^^(n) ^ 1, the detector network has very poor ability to disentangle the two polarization 
components of the incident waves, both of which are needed in order to determine D. As 
shown in Fig. ^, there are large regions on the sky in which the polarization sensitivity I — Ed 
of the LIGO/VIRGO network is poor, which correspond to the regions of high AD/D in 
Figs. and I fSll. 

Part of the reason for the low values of 1 — for the LIGO/VIRGO network is that the 
two LIGO detectors are nearly parallel, so that they access essentially a single polarization 
component of the gravitational wave field. [They were chosen in this way in order to enhance 
the reliability of detection of burst sources] . The addition of a fourth detector would greatly 
improve the polarization sensitivity of the network. In Fig. |^ we plot that fraction Q{e£,)/47r 
of the sky in which the polarization sensitivity is < 1 — Ed, for the LIGO/VIRGO network. 
We also as an illustration plot the same quantity for a hypothetical 4-detector network 
consisting of the LIGO and VIRGO detectors together with a detector in Perth, Australia, 
whose parameters [cf. Eq. (|90|) above] are assumed to be {9,ip,a) = (121°, 116°, 90°). 

The quantity T(n, v, ip) shown in Figs. |^ and ^ gives distance-measurement accuracy as a 
function of Earth-fixed coordinates {6, (p). The distance-measurement accuracy for coalescing 
binaries at a given right ascension and declination, averaged over many sources with different 
arrival times, will clearly be given by the average over of T (due to the Earth's rotation). 
Values of this averaged accuracy in the band 0.15 ^ | cos^^| ^ 0.65 are typically a factor 
of ~ 2 better than those outside this band, over the poles and near the celestial equator. 
Similarly, the average over ip of l/cr/)(0, p) is roughly proportional to the average maximum 
distance to which sources can be seen at a given declination; it does not vary by more than 
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~ 20 %. Note that the distribution of sources on the sky is expected to be approximately 
isotropic because the large distance 200 Mpc) to typical coalescences. 

Lower bounds for /\D / D can be obtained by combining Eqs. (0), (p6|), and (0) and 
minimizing over n, f , and ip. If we define 



o-max = max o-Dl^nj 



1.04 (for LIGO/VIRGO), 



(100) 
(101) 



we obtain the following lower bounds on AD/D: 



AD ^ 2D 

D ~ ^/o\^ To ' 
AD 1 



together with the upper bound for the overall SNR p 



(102) 
(103) 



D' 



(104) 



These bounds remain roughly valid when effects that are nonlinear in ro/-D are approximately 
taken into account [cf. Figs. ^^-[T^ below]. 



D. Extension of analysis to beyond 
the Gaussian approximation 

As explained in Appendix A, the Fisher matrix approach to calculating the probability 
distribution function (PDF) for the measured parameters is an approximation whose validity 
depends in part on the particular set of variables one uses to evaluate the Fisher matrix. In 
particular, the approximation works best for parameters on which the signal h(t) depends 
linearly. The key idea for dealing with the degeneracy limit f ^ 1 is to calculate the 
Gaussian probability distribution for the amplitudes Aa, which is exact because the signal 
depends linearly on these amplitudes (see Sec. [A6| below). Substituting Eq. ( [7I| ) into this 
PDF then yields the exact, non-Gaussian distribution for the parameters D, v, ip and (pc 
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where we mean "exact" in the context of the Markovic approximation An = 0. From this 
non-Gaussian distribution, values of AD can be determined which are more accurate than 
those given by Eq. (p6| ) in the regime v 1 and for low signal-to-noise ratios. 

In this subsection we calculate so-called Bayesian errors instead of frequentist errors. The 
distinction is carefully explained in Sec. [A 2| of Appendix A. The distinction is important only 
beyond leading order in 1/p, and hence unimportant elsewhere in this paper. In practical 
terms, the use of Bayesian errors means that the rms errors will be expressed as functions 
of the measured, best-fit values for the source parameters, instead of their true values. 

By using Eqs. ([7^ ) and (182|) one finds that the exponential factor in Eq. (|A10|) , given a 



gravitational wave measurement, is proportional to 

(105) 



exp 



Here the quantities Aa are the amplitudes that we measure at each detector (by using 
matched filtering). The corresponding PDF for the intrinsic amplitudes Aa is, from Eqs. (|72D 
and (|ATO|) , 



-| (A-i^)*(A-iB)e^^ 



p{Aa) = ATp^^'HAa) (106) 
X exp 

where 

Aa^{@-')abF^A,k''\ (107) 

Here p^^\Aa) is our a priori PDF for the amplitude parameters, and A/" is a normalization 
constant. 

As an aside, Eq. (|107| ) provides us with the maximum-likelihood estimator D (in the 
constant n approximation) of the distance to the binary in terms of the measured amplitudes 
Aa- This is because Eq. (ffl]), re-expressed in terms of hatted quantities, may be inverted to 
determine Z) in terms of the ^^'s: 



D = , 108 
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where [cf. Eq. (^T]) above] 

1 r 



6? 



Ar + 1^x1- |^+ + ^xl (109) 



and 

/3' = ^[ii+r+iixr+iiu^xi]- (110) 

The PDF p^^\Aa) in Eq. (|0|) represents our a priori information about the distribution 
of the parameters Aa [or equivalently from Eq. (|7T| ) of the parameters (D, f , (pc) ], given 
that a signal has been detected. Since we expect sources to be uniformly distributed in 
orientation and in space (on the relevant scales of ^ 100 Mpc), we take 

xe{D)e{D^,,-D)D^dD. (Ill) 

Here 9 is the step function, and the cutoff for distances greater than -Dmax is a (somewhat 
crude) representation of our knowledge that very distant sources would not have been de- 
tected. A suitable choice for -Dmax is the distance tq, cf. Eqs. (^Tj) and ( PB| ) above. Our results 
below are insensitive to the exact location of this cutoff, but it must be included to make the 
PDF (|106|) formally normalizable. Now let Dq, vq, ipo and 0co be the parameters obtained 



from the amplitudes A a by inverting Eq. (|7TD , so that, in particular, Dq = D. Substituting 
Eqs. ([71|) and ( [1 1 1| ) into ( |106|) yields a non-Gaussian PDF for the variables {D, v, ip, 0c) 
which depends on the parameters {Dq, vq, ipo, 4>co)- From this PDF it is straightforward in 
principle to calculate AD, by first integrating over v, ip, and 0c to determine the reduced 
PDF p{D) for D alone. If one first expands the argument of the exponential to second order 
in the quantities D — Dq, v — Vq, ip — ipo, and 0c — 0co, the result obtained is just Eq. ( pUf ) 
above, which is accurate to linear order in 1/p. 

Thus, in order to go beyond this linear approximation, one has to integrate the PDF 
( |106| ) over v, ip and 0c. Because this is difficult to do exactly, we now make an approximation 
which treats the correlations between D and {ip^ip^ to linear order in 1/p, but treats more 
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precisely the strong correlations between D and v. This approximation should give rough 
estimates of effects that are nonlinear in 1/p, and moreover removes the singularity in our 
previous result ( p6| ) at f = 1. The approximation consists of expanding the argument of 

- (j)co, and integrating over ip 



the exponential in Eq. ( |106| ) to second order in tjj — ipo and (f)c 
and 0c- One then obtains a function of v, D, vq, Dq, and ipo', the dependence on 0co drops 
out. This function is of the form (prefactor) x (exponential factor). The prefactor depends 
only weakly on D and v in comparison to the exponential factor, so we can approximate 



it to be constant [£2[- We then obtain the following PDF, which may also be obtained 



by substituting the transformation (0) into the Gaussian PDF for the variables a, (3 that 
corresponds to the variance-covariance matrix (0). 

The result is, in terms of the rescaled distance V = D/Dq, 

2 



dpiv, V) = MV^ exp I - ^ (^^ - 



2Ai 



n 2 



2V 2 

xQ{V)Q{D^^J Do - V)Q{1 - v^)dvdV. 

Here B is the step function. A/" is a normalization constant. 



(112) 



1 + Ed cos(4'0o) 



ro \ anil - e\)) 



I- Ed cos(4?/'o) 



(113a) 
(113b) 



and "^o = ^0 + A'?/'(n) [cf. Eq. ( PU[ ) above]. In terms of these quantities, the previous, 
approximate result (^Bj) is 



AD '2jvlAl + Ai 



l-vl 



(114) 



From the PDF ( |112|) one can numerically calculate the reduced PDF for T> alone. 



p{V) 



dvp{v, V), 



(115) 



54 



and thereby determine /S.D. As an example we show in Fig. ^ a plot of p{D) for a particular 
choice of the parameters Dq, vq, i/jq, and for a particular direction on the sky. The non- 
Gaussian fall off in this figure at large values of D is a general feature, although its magnitude 
in this example is larger than is typical. It can be seen that the distance measurement 
accuracy is a factor of ~ 2 worse than that predicted by Eq. (PB). 



Now from Fig. [1^ it can be seen that the value of D which maximizes p{D) is not the 
same as Dq, i.e., the Z)-component of the point {vq, -Do) which maximizes p{v, D). Hence, the 
"maximum-likelihood" method for estimating signal parameters is ambiguous — the results 
obtained for one variable depend on whether or not other variables are integrated out before 
the maximum is taken. As explained in Appendix A, we advocate as the "best-fit" value of 
D the expected value 

{D) = J Dp{D)dD = J Dp{v,D)dvdD, (116) 

instead of the maximum-likelihood estimate Dq. [Maximum-likelihood estimation will need 
to be used, however, to obtain initial estimates of the signal parameters]. Correspondingly, 
to estimate distance measurement errors we use the quantity 



AD ^ - ip? 

D ~ {V) ' ^ ' 

This can be calculated numerically from Eqs. ( |112| ) and (|115|) , and in general will depend in 



a complicated way on the parameters Ai, A2, and Vq, and very weakly on the rescaled cutoff 
Draax/Do- FoT the binary merger example of Fig. |10|, we show in Fig. |lT] how the accuracy 
( |117|) varies with vq, and in Fig. p!2| how it varies (through the parameters Ai and A2) with 
the distance Dq. 

The merger of a BH-NS binary of masses IOMq and 1.4Mq would produce a signal whose 
amplitude is 2.11 times stronger than the NS-NS merger of Fig. |12| [from Eq. ( ^3]) above]. 
Hence, taking also into account a cosmological enhancement factor of (1 + z)^^^ B, a plot 



of AD I D versus Dq for a BH-NS binary otherwise the same as the binary in Fig. [1^ would 
look roughly the same as Fig. |T^, but rescaled to extend to luminosity distances ~ 2 Gpc 
(the exact value depending on the cosmological model) 
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E. Simulation of what LIGO/VIRGO will measure 



In order to explore more completely the distance measurement accuracy ( 117 ) over 
the whole parameter space, we carried out the following Monte-Carlo calculation. Ran- 
dom values of Dq, vq, ipo, 0, and were chosen, distributed according to the measure 
(IDq dvQ dipQ d cos 6 dip . Those parameter choices for which the combined SNR ( [5^ ) (for NS- 
NS binaries) was less than the threshold of 8.5 were discarded, and samples were gener- 
ated until 1000 NS-NS signals had been "detected." Because of this thresholding procedure 
(which roughly corresponds to the actual thresolding procedure that will be used), the distri- 
bution of values of Dq, vq etc. for detected events will not be given by dDQ dvo dipo d cos 9 dip. 
For example, there is a significant bias in detected events towards high values of Vq, i.e., 
towards face-on binaries. 

Scatter plots of the distances Dq, signal-to-noise ratios p, and distance measurement 
accuracies (|117|) for these randomly generated data points are shown in Figs. |1^ - ^ We 
used the LIGO/VIRGO network functions shown in Figs. |^ and ^. Figures |13| - |1^ give 
some idea of the potential capability of the LIGO/VIRGO network. The distance scale in 
these graphs is determined by the detector sensitivity level (^ that we have assumed, which 
is uncertain to within a factor of ~ 2. The distance scale would also be 2 to 3 times larger 
for NS-BH binaries, as mentioned above. By contrast, the distribution of measurement 



accuracies, which we show in Fig. |Ty, is independent of the scale of the detector noise. This 
figure shows that the measurement accuracy will be better than 30% for over half of the 
detected sources. 

A relatively large fraction, about 1/5, of detected events have poor (> 50%) measurement 
accuracies. This is primarily due to the effect discussed in Sec. |IV C| : low values of the 
detector network polarization sensitivity 1 — ^^(n) over much of the sky. The effect of the 
polarization sensitivity can be clearly seen in Fig. |1^ which is a scatter plot of polarization 
sensitivity versus distance measurement accuracy. 

Finally, we emphasize that our results should be regarded as fairly rough estimates. 
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because we have neglected the following effects: (i) The spin-related modulation of the am- 
plitudes Aa mentioned in in Sec. |1V A| and discussed in Ref. (ii) the correlations between 



the variables D,v and ^,(f)c, except to linear order in 1/p [Q; and (iii) the correlations be- 
tween the parameters D,v,ip,(f)c and the "phase parameters" M., (3. As discussed above, 
we show in Appendix C that these correlations vanish to linear order in but there will 
be some correlation effects at higher order. Despite these neglected effects, we feel that the 
approximation method that we have used based on Eqs. (|112|) and (|115|) gives results that 



are more considerably more accurate than previous linear treatments [as summarized by 
Eq. (|96|)], because the dominant correlations at linear order in 1/p are those between D and 
f , and we have treated these correlations exactly. 



V. CONCLUSIONS 



Modulo the caveats in Sec. [Ill C|, we have confirmed the general conclusion that one can 



measure the binary's chirp mass M. with rather astonishing accuracy. While our estimates 
of /\M./M. are a factor of ~ 20 greater than those obtained from the less accurate New- 
tonian analysis [Il4| , |l5i , we have found that /\M./ M. should still be 0.01% — 1% for typical 



measurements. 

We have investigated the idea that detailed phase information might also allow accurate 
determination of the binary's reduced mass p. A calculation that neglected the effects of 
the bodies' spins on the waveform suggested that p might typically be measured to within 
~ 1%. However a more complete analysis showed that errors in p can be substantially 
masked by compensating errors in the spin parameter (3. Including the correlations with /3, 
we estimated that Ap/p 10% for low-mass (NS-NS) binaries and that Ap/p 50% for 
high-mass (BH-BH) binaries. Moreover, AMi/Mi and AM2/M2 are generally much greater 
than Ap/p unless M1/M2 > 1 (BH-NS case). 

These results are somewhat disappointing; it would have been more exciting to find that 
post-Newtonian effects allow both masses to be determined to within a few percent. In this 
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regard, however, it is useful to keep two points in mind. First, since typical measurements 
will have S/N ^ 10, one should detect events with S/N > 50 roughly 1% of the time. For the 
advanced-detector noise curve (H), and assuming the coalescence rates estimated in Ref. 0, 
such strong events should be seen ~ once per year for NS-NS binaries, and several times per 
year for NS-BH and BH-BH binaries [22|. For these strongest sources, measurement errors 



will be a factor of ~ 5 lower than their typical values. Second, the measurement-derived 
PDF on the parameter space constrains the values of Mi, M2, and /3 much more strongly 
than is indicated by their individual variances, as illustrated in Fig. ^ above. The large rms 
errors are due to correlations between the measured parameters; certain linear combinations 
of the parameters (eigenvectors of the covariance matrix) are determined with high accuracy 
40| . This may be useful when combined with information obtained by other means. 

With regard to potential accuracy of distance measurements, our key conclusions are the 
following: 

(i) We have confirmed the general conclusion reached previously [p!5|,p| that correlations 



between the distance D and other angular variables (primarily the angle of inclination of 
the binaries orbit) will reduce AD by a factor of typically 2 or 3 from the naive estimate 
AD/D = 1 /(signal-to- noise); see Fig. above. 

(ii) Distance measurement accuracy will depend strongly on the direction towards the 
source relative to the detectors, as shown in Figs. ^ and §. This is because of the different 
polarization sensitivities of the detector network in different directions, and the fact that the 
complex amplitudes of both polarization components of the incident waves are needed in 
order to determine the distance. The polarization sensitivity of the LIGO/VIRGO network 
is somewhat poor in this regard (because the two LIGO detectors are almost parallel); 
it would be substantially improved by the addition of a fourth detector. This provides 
additional motivation for the construction of additional interferometers around the world, 
which would also improve the angular accuracy of sky-location measurements 

(iii) Previous estimates of distance measurement errors have been accurate only to linear 
order in 1/D. Our results indicate that this linear approximation will be inadequate for 
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typical detected signals, so that the incorporation of non-linear effects will be necessary 
in order to accurately ascertain measurement errors [and also to accurately estimate the 
distances themselves; see Sec. [A3| of Appendix A] . 



(iv) We have carried out a Monte-Carlo simulation of distance measurement accuracies 
for a large number of randomly chosen sources, using a method of calculation which roughly 
estimates the non-linear effects, and incorporating the amplitude sensitivity and polariza- 
tion sensitivity of the LIGO/VIRGO detector network. Our results suggest that ~ 8% of 
measured distances will be accurate to better than ~ 15%, and that ~ 60% of them will be 
accurate to better than 30%. 
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APPENDIX A: ESTIMATION 
OF SIGNAL PARAMETERS 



In this appendix we review some aspects of the statistical theory of estimation of signal 
parameters as applied to gravitational wave astronomy. This subject has been concisely 
summarized in Appendix A of Ref . , and has recently been treated in detail by Finn |^ . 
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Hence in many places we merely write down the key results, without proof, in order to es- 
tablish notation and equations for use in the text. However, we also present some extensions 



to the formalism developed by Finn We carefully distinguish between Bayesian and 
frequentist estimates of errors, and discuss the validity of these two methods of error cal- 
culation. We show that maximum-likelihood parameter estimation, while useful, is not the 
optimal data-processing strategy, and, following Davis suggest the use of the so-called 
Bayes estimator. We derive an expression for the minimum signal-to- noise ratio (S'/A^)min 
necessary in order that the usual Gaussian approximation for estimation of measurement 
accuracy be valid, and explain how to treat degenerate points in parameter space at which 
the Gaussian approximation breaks down. Finally we give a discussion of the effects of 



including a priori information, which corrects the corresponding material in Ref. |26 



1. Basic Formulae 

The output of a network of detectors can be represented as a vector s(t) = 
{si{t), . . . , Sn^it)), where is the number of detectors, and Sa(t) is the strain amplitude 
read out from the ath detector. There will be two contributions to the detector output s{t) 
— the intrinsic detector noise n{t) (a vector random process), and the true gravitational 
wave signal h(t) (if present): 

s(t) = h(t)+n(t). (Al) 

We assume that the signal is a burst of known form, but depending on several unknown 
parameters 9 = {6^ , ... ,6^), so that h(t) = h(t; 0). Thus, we do not consider the cases of 
periodic or stochastic waves |12|. We also assume for simplicity that the detector noise is 
stationary and Gaussian. For the LIGO and VIRGO detectors, the stationarity assumption 
is justified for the analysis of short, burst waves |^. However, the actual noise may have 
important non-Gaussian components, the implications of which for the purposes of signal 
detection thresholds and data analysis are not yet fully understood. We do not deal with 
this issue here. 
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With these assumptions, the statistical properties of the detector noises can be described 
by the auto-correlation matrix 

CniT)ab = {Sait + r)Sfe(t)) - (^.(t + t)) 

= Mt + r)nfe(t)) - {uait + r)) (nf,(t)), (A2) 

where the angular brackets mean an ensemble average or a time average. The Fourier 
transform of the correlation matrix, multiplied by two, is the power spectral density matrix: 

/oo 
dre'-'f^CMab. (A3) 
-oo 

This satisfies the formal equation 

(naif) nbif'T) = l^if - nSn{f)ab, (A4) 
or more generally and precisely 

(expl^i J dtw{t) ■n(t)|^ =exp|-^ j^"c?/w^S„-w|, (A5) 

for any sufficiently well-behaved test functions Wa{t). Here tildes denote Fourier transforms, 
according to the convention that 

h{f) = / e^^'f'h{t)dt. 

J —oo 

We note that there are two different commonly used definitions of power spectral density 
in the literature. The above convention is used in Refs. |]3| JT^ , |26|j35| , |4^ . The alternative 
convention is to use a spectral noise density defined by S^p{f) = S'„(/)/2, as used in 



Refs. 121.27.56.57.58 



The Gaussian random process n(t) determines a natural inner product (. . . | . . .) and 
associated distance or norm on the space of functions h{t). As discussed in Sec. [II B| , this is 
defined so that the probability that the noise takes a specific value no(t) is 



p[n = no] oc e-("«l"°)/^ (A6) 
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and it is given by ||28 



-j^ ab y 



(g|h)^43fJ/ df UfT ^nif)-' hif), (AT) 

JO L J 

where 3? means "the real part of" . It also satisfies the equation 



((n|g) (n|h)) = (g|h), (A8) 

for any functions g and h. 

In this paper we are interested only in the estimation of signal parameters once a grav- 
itational wave burst has been detected. Thus, we suppose that we have measured some 
detector output s(t), and that it satisfies the appropriate criterion for us to conclude that 
it contains a signal of the form h{t; 9) for some unknown set of parameters 9: 

s(t) = h(t;^) +n(t). (A9) 

The central quantity of interest is then the probability distribution function (PDF) for 9 
given the output s(t). As Finn has shown pBf, this is given by 



p[9\s, detection] = M {9) 5 1 ^ W-s) _ (^10) 

Here M = A/'(s) is a normalization constant, and p^^\9) is the PDF that represents our a 
priori knowledge. 



2. Two types of measurement accuracy 

We now discuss how to characterize the accuracy of measurement of the parameters 
9. Normally statistical "one sigma" experimental errors are defined operationally in terms 
of the average of the actual errors over an ensemble of repeated identical measurements 
(which corresponds mathematically to the width of an appropriate PDF). Now in practice 
one cannot repeat or duplicate a given gravitational wave measurement, but in principle one 
can do so by waiting a sufficiently long time and throwing away all detected signals that 
do not match the original one. In this manner one can operationally define an ensemble of 

62 



"identical" measurements. The notion of error which results then depends in a crucial way 
on what is meant by "identical." One can either demand that the signals h(t) incident on the 
detectors be identical and consider the resulting spread in the values of the detector outputs 
s{t) given by Eq. ([X9|), or demand that the detector outputs s{t) be identical and consider 
the resulting spread in the values of the incident signals h(t). The two notions of error which 
result can be called Bayesian errors and frequentist errors, adopting the terminology from 
common usage in a more general context We now discuss in more detail the definition 
and meaning of these two types of error, in order to clarify the relationship between our 



method of calculating measurement error and previous work in this area [P6| , p7| , p!5| , p3| . The 
following discussion is based on that of Loredo . 

In the frequentist approach, one first specifies the algorithm the experimenters should 
use to determine the "best-fit" values 9 of the parameters 9 from the gravitational wave 
measurement s: 

^ = ^(s). (All) 



This is also called a statistic or estimator. Next, one assumes that Eq. ( |A9| ) holds for some 
value of 9, and by substituting this equation into Eq. (|A11|) , and using Eq. (|A6|) , one derives 



the PDF p{9 \ 9) for 9 given 9. Then the expected value with respect to this PDF of 9'' — 9\ 

U = 0') - 9\ 

9'p{9\9) d9 - 9\ (A12) 

gives the "bias" 6* of the estimator ^(s). The diagonal elements of the expected value of 
— 6'*) {9^ — 9^) characterize the measurement error. More specifically, we define 

S^REQ = Si^REQ[^~;^^(-)] (A13) 

'{^*[h(^) + n] - 9'] {9^[\i{9) + n] - 9^]' 

Here the notation on the first line indicates that Sfreq depends on the functional form of 
the estimator 9 as well as the assumed signal parameters 9, and the angular brackets on the 
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second line denote expectation value with respect to the noise n. The matrix ( [A13| ) is a 
measure of parameter-extraction accuracy that includes the effect of the bias, since 

^'^nEQ = {S0'S9^) + b^V, (A14) 

where 66' = - {§'). 

The physical meaning of the quantity (|A13|) is the following. Suppose that a large 



number of identical gravitational wave trains, described by the parameter values 6, impinge 
on the detector network. For each measured signal, the experimenters calculate using the 
algorithm 9 the best-fit values of the source parameters. Then the rms average deviation 
of these best-fit values from the true value 6 is given by Eq. ( |A13| ). Moreover, the usual 



method of implementing a Monte-Carlo simulation of the measurement process would also 
predict errors given by ( |A13D p3[ ]. 



In the frequentist method, one focuses attention on a particular incident signal h{t]6), 
and considers different possible measured detector outputs s(t). By contrast, in the Bayesian 
approach, one focuses attention on a particular measured detector output s. The error in 
measurement is simply taken to be the width (variance-covariance matrix) of the PDF (|A10|) 
for the true value 6 of 6 given the measurement s. Thus, 

^BAYES = ^BAYEs[^!P''°''(')] (^15) 



where (^*) = / 6' p{6 \ s)d6. Note that this measure of error depends on different quantities 
than its frequentist counterpart ( |A13|) — the measured signal s, and the a priori PDF p^^^ . 
The physical meaning of the quantity (|A15|) is the following. Suppose that a large number 



of different gravitational wave trains are incident upon the detector network, where the 
distribution of the wave parameters 9 is given by the PDF p^^^ . Only a small fraction of these 
will produce, at the output of the detectors, the signal s(t). In this small fraction, however, 
there will be some spread of values of the parameters 6, because of different realizations of 
the detector noise n{t) that combine with the incident waves to produce the measured signal 
according to Eq. (|A9|) . This spread is characterized by the matrix (|A15|) . 
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The measure of error ([A15|) characterizes the total amount of information that is con- 



tained in the measured signal s, which is independent of how the experimenters choose 
to process this signal. In practical situations, however, one typically would like to know 
what accuracy can be achieved by a given, imperfect, data-processing algorithm (e.g., one 
which takes a manageable amount of computer time). It is possible to define a more general 
Bayesian error that is appropriate for the situation where a particular algorithm or statistic 
6{-) is chosen to estimate the signal parameters 6 from the measured signal s. This measure 
of error is 

SLAVES = Si^AYEs[s,^(s);P^°H-)] (Al6) 

r - 0{s)) [6^ - e{s)) pie \ s)de. 



Physically this quantity is just the (square of the) rms average, over the small fraction of 
incident waves discussed above, of the difference between the true value 9 of the parameters 
and the "measured value" ^(s). It is clear that the rms errors Sg^YEsl^) ^(^)! "^ill be 

minimized and take on their minimum values Sg^yEsl^iP'^'^H')] when one chooses for 9 the 
so-called Bayes estimator 



9i,j,is)= 9'p{9\s)d9. (A17) 



One final point about Bayesian errors is the following. Suppose that the experimenters 
calculate from the measured signal s the best-fit value 9 = 9{s), and then discard all the 
remaining information contained in the signal s. Then there are very many signals s' that 
could have been measured and that are compatible with the experimenters' measurements, 
in the sense that ^(s') = 9 [Q. Correspondingly, there is a larger spread of possible values 



of 9, and hence the predicted rms measurement errors based on the measurement 9 alone 
are given by the following modification of Eq. ( |A16|) : 



^BAYES — ^BAYEs[^!P''°''(')] (^18) 

[9' -9) [9^ -9)p{9\9)d9. 
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Here p{6 \ 9) is the probability distribution introduced by Finn for the true parameter 



values 9 given the estimated values 9. It is given by the standard Bayesian formula 

p{9\9)=J\fp^^\9)p{9\9), (A19) 

where M = M{9) is a normalization constant that depends on 9. Note that the matrix ( |A18| ) 
depends only on the measured value 9 of the estimator and not on its functional form 9{-). 

The predicted measurement error (|A18| ) differs from the previously defined measurement 
error (|A16|) because the measured signal s contains information about the likely size of the 
error, so that discarding s makes a difference. For example, suppose that a detector-output 
data train contains a signal from a coalescing binary, and that by some standard algorithm 
the experimenters determine best-fit values of the binaries parameters. Then given these 
best-fit values, one can estimate the likely size of the measurement error — this is given by 
Eq. ( |A18|) . However, if they also determined that the data train contains an uncommonly 



large non-Gaussian burst of noise that accounts for 20% of the estimated signal amplitude, 
the estimates of the likely parameter-extraction errors would clearly have to be modified. 

Which of the above-defined measurement errors is appropriate to assess the capability 
of the LIGO/VIRGO detector network? It is generally accepted that, if one has a given 
measurement s, the Bayesian approach is the fundamental and correct one, and that the 
frequentist approach is justified only to the extent that it reproduces the results of Bayesian 
analyses. This is essentially because, given a particular measurement s, it is irrelevant to 



consider an ensemble of other, different measurements s' |Q. However, for our purpose 
of trying to anticipate the capability of gravitational wave detectors before any measure- 
ments are available, it seems that this message looses its bite. It certainly seems reasonable 
to imagine a fixed gravitational wave-train incident upon the detector network, and to in- 
quire about the spread ( [A 131 ) in measured values of the source parameters due to differing 



realizations of the detector noise. 

In fact, there is a certain sense in which frequentist errors and Bayesian errors are equiv- 
alent, which is well known: the average of the predicted frequentist error over the whole 

66 



parameter space is the same as a suitable average of the predicted Bayesian error. Thus, 
in a sense the same errors are being calculated in each case; it is just their dependence on 
parameters that is being changed. In particular, if the predicted errors do not vary strongly 
with the parameters 6', then the two types of error will be approximately equal. A precise 
statement of this "equality of averages" which is straightforward to derive is 

j^^^'iAYEshm;P^'\-)] (A20a) 
r rift 

= J j^^^tAYEs[0;p^'\-)]- (A20b) 

Here the various matrices S*-' are defined in Eqs. ( |A13| ), ( |A16| ), and ( [A18| ), respectively, and 



the factors Af{s) and J^{6) are the normalization constants appearing in Eqs. ( |A10| ) and 
( |A19| ), respectively. The (formal) measure Vs is defined such that 



(F[n]) = / F[n] e-Hn)/2^ (^21) 

for any functional F[n] of the noise n. 

Because of this equality of averages, we conclude that either Bayesian or frequentist errors 
can be used to anticipate the capabilities of LIGO/VIRGO, essentially because one is only 
interested in the range of possible errors and not their value at a fixed point in parameter 
space. Similarly, if one is using Bayesian errors, it is appropriate to use the matrix (|A18|) 



instead of ( [A16| ) to anticipate measurement accuracies, since from Eqs. (|A20|) the measure 
of error (|A18|) is simply an average of ([A16| ) over values of s for which ^(s) = ft. This 



conclusion has already been reached in a recent paper of Finn's p6[ in which he advocates 
the use of what in our notation is essentially Sg^yEsl^MLiP^^H')]' where 6'ml is the so-called 
maximum-likelihood estimator (see below). [However, his calculation of this quantity does 
not incorporate the a priori PDF quite correctly, as we show below.] Previous analyses of 
parameter-extraction accuracy for gravitational wave detectors by Echeverria |53] and by 



Krolak and collaborators P^ , |27|P5| , |55|] have used the frequentist error Epj^^ql^^; ^^ml(-)]- By 



contrast, in Sec. |I^ of this paper we have calculated the Bayesian error 
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Si|AYEs[s;P^°n-)] = SjfAYEs[s,^BE(s);p(°)(-)], (A22) 

because, as we argue below, it is more accurate to use ^be rather that ^ml- 

One final important point about the two types of error is the following well-known fact: 
to leading order in where p is the signal-to- noise ratio [Eq. ( |A27] ) below], the two 
approaches yield identical results. More specifically, assuming the Gaussian noise statistics 
we have 



S^req[^~; ^ml(-)] = sLyes[^ml;p(°H-)] X [l + O(p-i)] (A23) 

when ^ML = 6'. Moreover the same quantity [given by Eq. ( |A35| ) below] is also obtained 
to leading order using the Bayes estimator ( |A17| ), and also from the estimator-independent 



measure of error ( |A15|) . This is essentially because to this order, all the PDFs are Gaussian. 



These assertions are straightforward to prove using the tools developed by Finn |]26|, and 
moreover are well-known in more general statistical contexts. Hence the distinctions that we 
have been drawing only matter when effects that are non-linear in 1/p contribute significantly 
to the predicted accuracies (as for example when measuring distances to coalescing binaries), 
or when p is sufficiently small that the a priori information represented by p'-^^ becomes 
significant. [However, this may be the rule rather than the exception for typical detected 



gravitational wave bursts; see Sec. IV D above. 



3. Choice of data-processing algorithm 9{-) 

Given a particular measurement s, the PDF ( [AlOp in principle contains all the informa- 
tion contained in the measurement about the source parameters 9. However in practice one 
often wants to focus on a small portion of this information, by calculating a "best estimate" 
value ^(s) together with estimates of the statistical errors. The choice of estimator 9{ ) is 
determined both by practical considerations, and by whatever criteria are adopted to judge 
"good" estimators; there is no unique choice. 
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One obvious criterion is to choose that estimator which minimizes the expected error in 
parameter extraction. However, one could choose to minimize either Bayesian or frequentist 
errors, and also the errors depend on the parameter values {9 or 9). If one minimizes the 
average over parameter space of the measurement error [as given by the common value of 
Eqs. ( [A2U| )], then the resulting best estimator is just the Bayes estimator ( |A17| ), which we 



have used in Sec. |1V D[ Its use for gravitational wave data analysis has been suggested by 
Davis |^5|. Unfortunately, as Davis indicates, calculation of the Bayes estimator is very 
computationally intensive, as it typically involves a multidimensional integral of a function 
whose evaluation at each point requires the numerical calculation of an inner product of the 



type (^). Our application of the Bayes estimator in Sec. |IVD| was an exception in this 



regard, because all the inner products could be evaluated analytically. It seems likely that 
the Bayes estimator will be used only after preliminary estimates of the signal parameters 
have been made using Wiener optimal filtering. The use of the Bayes estimator also goes 
by the name of "non-linear filtering" p5[ . 

A simpler estimator that has been proposed by Finn pB|JT^ , Krolak and others in 



the gravitational wave data-analysis context is the so-called maximum-likelihood estimator 
^ml(s). This defined to be the value of 9 which maximizes the PDF ( |A1(J| ). It is convenient 



because it is closely related to the Wiener optimal filtering method |jT2| that will be used 



to detect the signals — the detection procedure outlined in Sec. | will essentially return 
the maximum-likelihood estimates of the source parameters (see below). However, once the 
signals have been detected, there is no reason to only use maximum-likelihood estimation 
— other estimation methods can be used to give better results. Hence, the quantities 
^bayesI^ml;^^"^')] ^FREQ[^j ^ml(-)] represent the potential accuracy of measurements 
only to leading order in 1/p. (If maximum-likelihood estimation is the only estimation 
method used, then they represent the actual accuracy of measurement). We note that the 
quantities ^ml(s) and 9be{s) can differ by substantial factors for detected gravitational wave 
signals, as for example in Fig. above where (D) = IAADq. 

Maximum-likelihood estimation also has the following disadvantages. First, as discussed 
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in Sec. [IV D| , the maximum-likelihood estimator for a particular variable does not necessarily 
maximize the reduced PDF for that variable obtained by integrating over the other variables. 
By contrast, the value of the Bayes estimator ( [A17| ) for a given variable does not depend on 
whether or not other variables have been integrated out. Second, the best-fit point obtained 
by the maximum-likelihood method depends on the choice of variables used to parametrize 
the waveform h{t; 9). For example, in Sec. |T| of this paper we could have used as variables 
either the individual masses Mi and M2 of the binaries components, or the chirp and reduced 
masses M. and /i. Since probability distributions for (Mi,M2) and (Ai,^) are related by a 
non-trivial Jacobian factor, a local maximum of one of them will not correspond to a local 
maximum of the other. A slightly different kind of maximum-likelihood estimator, which 
maximizes the likelihood ratio A(^^) oc exp [ — (h(^^) — s | h(^^) — s)/2] instead of the PDF 
( |A1CI|) , does not suffer from this problem. This is the maximum-likelihood estimator that is 
usually discussed in the statistics literature. However, it does not take into account in any 
way our a priori knowledge. 

We conclude that calculations of measurement accuracy using 6'ml represent the true 
potential measurement accuracy only to leading order in 1/p. If this leading order approxi- 
mation becomes invalid (as occurs for sufficiently small SNR's), then one should use instead 
either the Bayesian error (|A22|) or ^bayes[^^EjP^^\')]- could also use Spj^p^q[6'; 6'be(-)]) 
but this is much more difficult to calculate than ([A22|) when the large p limit does not apply. 



4. Relation between maximum-likelihood estimation 
and Wiener optimal filtering 

In Sec. H, we discussed a method for finding best-fit parameters 6 which was based 
on maximizing the overlap of the measured signal with theoretical templates [cf. Eq. 
above]. We now briefiy indicate the relationship of this method to the maximum-likelihood 
procedure. That the two methods are equivalent in general has been shown by Echeverria 
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Given the measured signal s(t), define for any 6 tlie quantity 



'(h(C)|h(0)) 



Tliis is tlie signal-to- noise ratio (SNR) defined in Eq. @, and can be calculated by integrating 
the signal s(t) against a Wiener optimal filter whose Fourier transform is proportional to 
Sn(/)~^ ■h(/; 9). The quantity p[0\ is a random variable with Gaussian PDF of unit variance. 
Its expected value is zero if no signal is present, when s{t) = n(t). If a signal is present, so 
that Eq. (|A9|) holds for some 9, then the expected value of p[9] is 



Now if the a priori probability p^^\9) can be approximated to be constant, then the value 
^ML of 9 which maximizes the PDF ( [A10| ) for a given signal s also maximizes p[9] Hence 
we can find 9ml (up to the overall amplitude of the signal [^) by computing the overlap 
( |A24|) of the signal with various templates, and by choosing the template which gives the 
maximum overlap. 

When a signal is present and the signal-to- noise ratio is large, the maximum value p[^ml] 
of p[9] will approximately given by 

p[9ml? ~ (h(^~) I h{9)) ^ (h(^ML) I h(^ML)). (A26) 
The quantity 

p' = ihi9)\h{9)) (A27) 

is what is usually referred to as the (square of) the SNR of the signal h(t; 9). When correlated 
sources of noise are unimportant so that the matrix ( |A3| ) is diagonal, this overall SNR will 
be given by combining in quadrature the SNR's for each individual detector, cf. Eq. (H) 
above. 
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5. The Gaussian approximation 
and conditions for its validity 



We now consider the high signal-to-noise hmit in which many of the subtleties that 
we have been discussing become unimportant. In particular, in this limit the Bayes and 
maximum-likelihood estimators become identical. From Eq. ( AlU ), the maximum-likelihood 
estimator 9ml satisfies 

(h,(^ML) I h(^ML) - s) - [lnp(o)],(^ML) = 0, (A28) 

where the subscript , i means derivative with respect to for 1 < i < k, and k is the number 
of parameters. If the a priori information is unimportant so that the last term in Eq. ( |A28| ) is 



negligible, then as outlined in Sec. JTT| the following simple geometric interpretation applies: 
Let S be the finite dimensional surface formed by the set of all signals h{t; 9) in the space of 
all possible signals h()f:). Then the measured signal s(t) will generally not lie on the surface 
iS, and the best-fit point h(t; ^ml) is just that point on S that is closest to s(t), where 
distance is measured using the norm ||f|p = (f | f ) derived from the inner product ( |A7D . 
Correspondingly, h(6'ML) can be obtained by just dropping a perpendicular from s(t) onto 
the surface 5, which is the content of Eq. ([A28|) and is illustrated in Fig. |l]. 



When the SNR p is sufficiently large, one can find approximate expressions for 
Sfreq[6'; ^ml(')] Sbayes[^ml; Such a calculation has been carried out by Finn 

26| . We now briefiy outline the calculation, and also extend it to determine the next to 



leading order terms in an expansion in 1/p, in order to determine how large p needs to be for 
the leading order term to be a good approximation. Throughout this subsection we assume 
that the a priori PDF p'-^^ is approximately constant; in subsection |A 6| below we consider 
the effects of non-constant p^'^\ 



First we find an approximate solution to Eq. (|A28|) . Abbreviating 6'ml as 9, inserting 



Eq. ( [A9| ) into Eq. ( |A28|) and expanding in the difference 69 = 9 — 9, we obtain 



9' = 9' + 6^^^9' + 5^^^9' + 6^^^9' + O(n^). (A29) 
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Here 



S^'-^e^ = (m-'Y (n|h,), 



(A30) 



where 



(A31) 



is the so-called Fisher information matrix [cf. Eq. ( p!0D above]. The second-order term §^"^^9 



IS 



(A32) 



In this expression and below we have for brevity omitted the commas denoting derivatives, 
and all quantities are evaluated at 6. We lower and raise indices with the tensor ( |A31| ) and 
its inverse, so that, for example. 



1 N 



(A33) 



There is a similar but more complex expression for 5^^^ 9'^. 

Equations (|A29|) and ( [A^ ) now determine the PDF p{9\9). Using Eq. (|X8|) and its 
extension to fourth order moments, and Eqs. (|A29|) and ( |A13|) , we obtain 



(A34) 



At leading order, p{9\9) is a multivariate Gaussian with mean 9 and variance-covariance 
matrix proportional to given by the first term in Eq. ( |A34| ) : 



(A35) 



cf. Eq. above. The correction term oc in Eq. (|A34|) is 



= (h\ I h^'^) - (h\ I h,) (h^-'^lh') 
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I y\) - (K^ I h\) + (h^^' I hfc) (h'^ I h' , 



,fc I ul 



+ (h' I h^*=) (hfc I h',) + 2 (h^ I hn) (h'^ I h^') . 



(A36) 



In the case where there is only one variable so that 6 = {6^, . . . , 9^) = (0^), it follows from 
Eqs. ([A3l ) and ( |A36| ) that 



1 + 



15 (h" 


h')' (h' 


h'")" 


4 (h' 


h')^ (h' 





(A37) 



where primes denote derivatives with respect to 6^. The correction terms in the square 
brackets in this expression will be small whenever 



|h"|| < \\h 



'I |2 



(A38) 



and 



< llh'lP. 



(A39) 



Using Eq. (|A7|), the equation = (h | h), and assuming for simplicity that S„(/) = S'„(/)l, 
Eq. ( |A38| ) reduces to the condition 



^ (((h"t.h")/(ht.h))) 



(A40) 



(((h't-h')/(ht-h)))2' 

where for any function of frequency F{f), we define the weighted average {{F{f))) to be 
(Fh|h)/(h|h). 

Equations ([A38|) and ([A39| ) give sufficient conditions for the Gaussian approximation to 
be valid, when there is only one unknown parameter 6^. When there are several unknown 
parameters, a generalization of Eq. ( |A38| ) is obtained by interpreting the prime to mean the 
operator v^d/d9^ that differentiates in some direction in the space of parameters 9, and 
requiring the condition to hold for all directions f'. This yields the condition 



i\\2 



(A41) 
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which is required to hold for all w*. We note that, although Eq. ( |A34|) does correctly indicate 



the regime ( |A4CI| ) where the Gaussian approximation is valid, the correction term (^^S*-' is 
not an accurate expression for the leading order correction to the measurement accuracy, 
because as we have argued above the true potential measurement accuracy is given by using 
the estimator ^be(') and not ^ml(')- 

One frequent source of confusion about the leading order expression ( |A35|) for the mea- 
surement error is the following. A general theorem in statistics called the Cramer-Rao 
inequality [^,0 states that for any unbiased estimator 9 [|6^ , 



Sfreq[^;^(-)] > W) ■ (A42) 



Hence, one might imagine that the quantity ( |A35| ) is a lower bound for the accuracy ob 



tainable by most reasonable estimators, and also for low signal-to-noise ratios. That this is 
not the case can be seen from the following consideration, which we discuss in the body of 
the paper: at degenerate points 6q for which the signal derivatives dh/d6^ become linearly 
dependent, the matrix ( [A31D becomes degenerate, and the predicted rms errors given by the 



matrix ( |A35D become infinite. More careful calculations of, for example, Sfreq[6'; 6'ml(')] 
at such degenerate points, going beyond linear order, give finite results. Hence the inverse 
of the Fisher matrix is not a generic, useful lower bound. The reason that the Cramer-Rao 
inequality does not apply is that most estimators are not unbiased and cannot easily be 
made so. When one generalizes the inequality ([A42|) to incorporate the effects of bias ||57|| , 



an extra factor appears on the right-hand side multiplying the Fisher matrix, which can 
be small. This can allow Sfreq[6'; ^(O] be much smaller than the inverse of the Fisher 



matrix, for some statistics 6. 



6. Incorporation of a priori probabilities 

We now turn to the effects of a priori information. First, we remark that it is not 
necessary for a priori information to be very detailed or restrictive in order that it have a 
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significant effect on parameter-extraction accuracy. All that is necessary is that it be more 
restrictive than the information contained in the waveform, for some of the parameters 
6**. In other words it will be important whenever the statistical error ((A^^*)^) which we 
obtain from Eq. ( |A35| ) for some parameter 6'* is much larger than our a priori constraints 
on 6^. For example, this would be the case if we obtained rms errors for measurements of 
the dimensionless spin parameter a of a black hole to be larger than one, since we expect 
\a\ < 1 always. If we include such poorly determined variables in a calculation of the 
variance-covariance matrix S and neglect the a priori restrictions, then the results obtained 
for the rms error in ^* may be severely overestimated. This is not unexpected; what is more 
surprising is that due to the effects of correlations, the rms errors obtained for the other 
parameters may also be overestimated by large factors (see, e.g.. Sec. [Ill B| above). We 



now extend the approximate calculations of the previous subsection to incorporate a priori 
information, and also now calculate Bayesian as well as frequentist errors. Our results in 
this subsection correct Eq. (3.19) of Ref. [p6| . 

Roughly speaking, a priori information will be unimportant when the PDF does not 
vary substantially within one or two sigma of the best-fit point 6. This condition is logically 
independent of the condition ( [A4U| ) , although both will be satisfied in the high p limit. 
Hence, we can treat separately deviations from the leading order measurement accuracy 
( |A35|) that are due to second-order derivatives h of the signal [cf. Eq. |A36| above], and 



that are due to a priori information. In the remainder of this subsection we therefore assume 
the condition ( |A40|) and consistently neglect all second-order derivatives h jj. In particular 



we treat the Fisher matrix (|A31|) as a constant in this approximation. [Note that our results 



will be exact in the case where the dependence of h(0) on the parameters 6 is exactly linear, 
as in Sec. p.V D| above.] 

We start by considering the Bayes error ( |A15[ ). When we neglect second-order derivatives 



of h we find that the PDF ( [AlOp takes the form 



p[e I s, detection] = X' p^^\9) 
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X exp 



(A43) 



Here we have decomposed the measured signal according to 



s = s-^'hj + s"* 
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(A44) 



into the 



where {h.^i | s-*-^ = for 1 < i < /c, and have absorbed a factor of exp 
normahzation constant M' . If the PDF p^^^ simply restricts the allowed ranges of the param- 
eters, then the PDF ( |A43| ) is a truncated Gaussian distribution whose variance-covariance 
matrix Sbayes[s;p^°''(-)] will normally be within a factor of ~ 2 or so of T^. If p^^^ is 
approximately Gaussian with variance-covariance matrix So, then we see from Eqs. ( |A15| ) 
and flTHI) that 



SBAYEs[s;p(°)(-)] = {r + So^}"'. 



(A45) 



This is the formula which we use in Sec. [IIIB| above to incorporate our a priori knowledge 
about the spin parameter (3. 

Next we calculate an approximate expression for the second type of Bayesian error given 
by Eq. ( |A18|) , which is appropriate for the situation where we discard all information about 
the measured signal s except the best estimate values ^(s) of the parameters. For simplicity, 
we assume that p^^\6) is a Gaussian with mean 9q and width Sq. We also use the maximum- 
likelihood estimator 6'ml; however, the same results are obtained for the Bayes estimator ^be- 
From Eqs. (|A28D and (|A9|) and neglecting second-order derivatives of h, we find 



[K I n) + fSn ^ 



^0 i^Aol-e^), 



(A46) 



where = F + Sg . Together with Eq. (|A8|) this implies that 



p{9Mh\0) oc exp 



-v^ ■ F ■ V 



(A47) 



where 



v = ^-Si-r-^-Si-Soi-^o 

= Si ■ r ■ - const^ 
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(A48) 
(A49) 



Using Eqs. ( |A13| ), ( |A47| ) , and ( |A48|) we see that the result ( |A34]) becomes modified to read 



Sfreq[^;^ml(-)] = r"' + b®b, (A50) 

where the bias b = Si ■ Sg ^ ■ {9q — 9). A more interesting quantity is the Bayesian error 
( |A18|) , which from Eqs. ( |A19|) and ( |A49|) is given by 



SBAYEs[^ML;p(°n-)]"' = ^BAYEsIObe; P'-'H-T' (A51) 

= Soi + r-Si-r-Si-r. 

This expression gives approximately the same results as Eq. ([A45| ) , the differences never 
being more than ~ 25%. The variances S** given by Eq. ( |A51| ) are always larger than those 
given by Eq. (|A45|) , as a result of our having thrown away all the information in s apart 



from ^(s). 

The result ( |A51| ) disagrees with a corresponding analysis of Finn [Eq. (3.19) of Ref. 



The reason for the disagreement is that Finn solves Eq. ( |A28D to obtain as a function 



of ^ML and n, and then invokes the PDF ([A 61) of the noise to find p{6\9ml)- This method 



of calculation [analogous to the method used for calculating p{9ml\9)] is invalid because it 
implicitly assumes that 

p[n = no I ^ml] = p[n = no], (A52) 

which is not the case. The fact that Eq. ( |A52| ) does not hold can be seen from the joint 
PDF for 9, 9ml and n, which is 

p[9ML,0,n] ocpW(^~)e-("l")/2 

x5(^ML-^ML[h(^) + n]). (A53) 



7. Treatment of degenerate variables 



As we have noted in Sec. [IV D| , the accuracy of the linear approximation ( [A41D which 
yields the simple PDF (|A43|) , depends in part on what set of variables 6'* are used in the 
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calculation. A different PDF will be obtained from this approximation if one first makes 
a non-linear change of co-ordinates 9^ 6^ (6^). Hence, the PDF ([A43|) will approximate 



most closely the true PDF when it is computed using variables for which h^ij is as small as 
possible. 

Consequently, there are two qualitatively different ways in which the linear approximation 
may break down. First, for sufficiently low signal-to-noise ratios, the extrinsic curvature of 
the surface S formed by the set of waveforms h{t; 6) may be sufficiently large that Eq. (|A41|) 



is not a good approximation for any set of coordinates 6^. In this case the "Gaussian" method 
breaks down completely. Second, the approximation may break down simply because of a 
bad choice of coordinates. This is usually the case at points of degeneracy where the vectors 



dh/dO^ become linearly dependent, which we discuss in Sees. |T| and |^ and at the end 
of Sec. [A 5| above. At such points the straightforward linear approximation method breaks 
down, but frequently one can find a non-linear nonlinear change of variables of the form 

6' = e\e^) (A54) 

such that the vectors dh/dd^ are not linearly dependent. One then obtains from from the 
linear approximation a Gaussian PDF in the variables 6^. Substituting the relation (|A54[) 



into this PDF and multiplying by the appropriate Jacobian factor gives a non-Gaussian PDF 
in terms of the variables 6'* (as in Sec. |1V D| above). From this PDF, measurement errors for 
the variables 6^ can be calculated. As has been pointed out by Markovic 0, measurement 
errors at degenerate points in parameter space typically scale like 1 / ^/p instead of like 1/ p. 
This is true if the lowest order derivative of h which is non- vanishing in all directions is the 
second derivative, as can be seen from, e.g., Eq. (|A29|) above. 



APPENDIX B: APPROXIMATE CONSTANCY OF 
THE SPIN PARAMETER (3 THAT INFLUENCES 
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THE WAVEFORM'S PHASE 



The leading order contribution of the bodies' spins to the secular growth of the 



gravitational-wave phase has been derived by Kidder, Will and Wiseman |iT], and is given 
by the term proportional to in — P in Eq. (^91) . The quantity /? is defined by Eq. (|50D and 
depends on the masses of the two bodies Mi and M2, their spins 5*1 and S2, and the unit 
vector in the direction of the orbital angular momentum L. Over the course of the inspiral /? 
will evolve, because the directions of the vectors L, Si, S2 will be changing due to spin-orbit 
and spin-spin interactions. 

Nevertheless, in our analysis in the body of the paper, we have assumed that the factor 
X = 4:71 — (3 which appears in Eq. ( ^9]) can be treated as constant. This assumption is 
necessary to make the analysis tractable. In this appendix, we present evidence which 
strongly suggests that x is always constant apart from some small amplitude oscillations, 
showing that our assumption of constant x is a reasonable one for all coalescing binaries. 
We calculate the evolution of x by integrating the orbit-averaged equations (0) governing 
the evolution of the spins, using both analytic and numerical methods. A more complete 
discussion of the evolution of the spins and orbital angular momentum, and of their influence 



on the emitted gravitational waves, can be found in Ref. ||^ . 

We start by introducing some dimensionless variables. Let Sj be the unit vector in the 
direction of Sj for j = 1,2, and define 

ai = Si- L, (Bl) 
«2 = ^2 ■ L, (B2) 
«3 = 5i X ^2 ■ L, (B3) 

and 

a4 = ^1 ■ ^2. (B4) 
The aj's are not all independent variables as they satisfy the constraint 
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(B5) 



The reason that it is convenient to use these variables is the following. The spin-evolution 
equations ( |51| ) comprise nine equations in nine unknowns, with three conserved quantities 
(the magnitudes of the three vectors). Thus, there are effectively six degrees of freedom. If 
we specify the three independent values of the variables ai, . . . , ^4, then the remaining three 
degrees of freedom can be parametrized by an overall rotation matrix. More precisely, given 
the vectors 5*1, S2 and L, there will be a unique rotation matrix R which takes L into L! = 
(the unit vector along the z axis), 5*1 into a vector 5*^ in the x-z plane, and 5*2 into some 
S2. The vectors 5*1, 5*2 and L are determined by R and by the variables ai, . . . , 04. Hence, 
the variables <yj{t) for 1 < j < 4 and R(t) can be used instead of the vectors themselves to 
parameterize a solution to the spin-evolution equations. Now it turns out that the evolution 
of the ctj-'s decouples from the evolution of R, in the sense that each daj/dt depends only 
on ai, . . . , ^4 and is independent of R. This greatly simplifies our analysis. 

If we use units in which M = 1 and define Sj = \Sj\ for j = 1,2, then we obtain from 
Eqs. ([5T|), ( |B1| ) and (p!3| ) the following coupled system of equations for ai, . . . , a^: 



da I 


-15 ■ 


dr 


~ 128/i [ 


da2 


15 ■ 


dr 


~ 128/i [ 


da4 


— 15q;3 


dr 


128/i 



1 SiOi 

M2 ^ 

1 S2Q;2 

M2 Ml 
Ml M2 



32(^3 

sias 

L + Siai — 52^2 



(B6) 
(B7) 
(B8) 



Here L = yU-^/r denotes the magnitude of the orbital angular momentum, and we have 
changed the dependent variable from time t to orbital separation r. The omitted equation for 
da^/dr can be obtained by combining Eqs. ( p5|) - (|B8|) . From Eq. (|50|), the spin parameter 
f3 is given in terms of these variables by 
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12 



[SiOi + S2tt2) 



+ 



25 



AM1M2 

We have numerically integrated the equations 
angular momentum directions, for the cases of NS-NS, NS-BH and BH-BH binaries. We 



(B9) 

(PSP for various initial spin and 



assumed all neutron stars have masses of 1.4 Mq, and black holes have masses 10 Mq. 
We integrated inwards, starting at that value of r at which the emitted waves enter the 
LIGO/VIRGO waveband at 10 Hz, and ending at r = 6M near the last stable circular orbit 
pO| . In the special case that one of the spins vanishes, it can be seen from Eqs. (|^) and 
(pT|) that f3 will be conserved [Q. Hence we took both spins to be non- vanishing. We also 



assumed that their magnitudes are maximal, so that sj = Mj for j = 1,2 [cf. Sec. pi B|] , as 
these are the values which can be expected to give the largest changes in f3. 

Typical results are shown in Figs. ^ - The factor x = 47r — /? undergoes small 
oscillations with an amplitude of order 0.1 which is small compared to the mean value of x- 
This mean value depends on the mass ratio and on the initial spin directions, but always 
lies between An — /5max ~ 3 and An + P^aax ~ 22, where Praax is as given in Sec. piB| . The 
angles between the vectors given by ai, . . . ,04 also oscillate, all with the same frequency. 
[This frequency is not the frequency with which the total spin S = Si + S2 precesses around 
L as that precession does not change the angles between the vectors, and thus is not 
described by Eqs. (|B§) - (IBg])]. 



Some insight into the behavior of the general solutions to Eqs. ( PB| ) - ( |B(j| ) can be gained 
by considering the special case when the magnitude of one of the spins (say Si) is small, so 
that Si = 15*1 l/M^ ^ 1. This condition will sometimes be satisfied by NS-NS and BH-BH 
binaries, but will always be satisfied by NS-BH binaries since all compact bodies satisfy 
\Sj\ ^ Mj. Below we find analytic solutions to first order in Si. As we now describe, in 
the approximation si <^ 1 the amplitude of the oscillations of f3 (and hence also of x) is 
always smaller than ~ 1/4, for all initial spin directions and for all mass-ratios. Although 
rapidly spinning NS-NS and BH-BH binaries will not satisfy si <^ 1, nevertheless we find 
that amplitudes of the oscillations of (3 in the numerical solutions agree roughly with those 
predicted by the small spin approximation. [For some special initial spin directions, such as 
ai = 0(3 = 0, the analytic solutions are poor approximations to the numerical solutions, but 
in all such cases that we have checked, the amplitudes of the (3 oscillations are still ^ 0.2] 

The solutions to first order in si can be written as 
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a^.(r) = af\r) + af\r) s, + 0{si), (BIO) 



for 1 < J < 4. Now as we have already mentioned, it can be seen from Eqs. (^OD and (|5TD 
that when si = 0, the angle between 5*2 and L is conserved, so that (3 is constant. However, 
in this case the angles between the small spin 5*1 and the other two vectors will not be 
conserved. Thus, the zeroth order solutions a^^^ will be non-constant. We start by deriving 
these solutions. 

Substituting si = into Eqs. ( p6|) - (p8| ) we find that 02°'' is constant, i.e., 0*2^ = 02,1 = 
0^2 (^i) 5 where rj is the initial orbital separation, and that 

^ - -haf^ (Bll) 



dr 



(0) 

/i4ar- (B12) 



(ia4 _ (0) 



dr 
Here 

h = (B13) 
128/iM2 ^ ^ 

15 

h, = ^{s,a2, + L6), (B14) 

and S = {Mf — M|)/(MiM2). The coefficient /14 is non-constant as L = fl^/r depends on 
r. However, since it will typically vary slowly compared to the oscillations in the angles, we 
can approximate it to be constant. [The evolution of gives rise to a slow evolution in the 
amplitude and frequency of the oscillations of the a;j°^'s]. Defining 

a± = h^a'^^ ± hiaf^ (B15) 

we find that a;_(r) is constant, a_(r) = a_^j = a„(rj), and 

^ = -2hh,aP. (B16) 

This equation can be solved by combining it with the constraint (p5|) . To zeroth order in 
Si, the constraint can be expressed using Eq. ( P15|) in the form 

af\ry + u^[a4r)-a+f = K^, (B17) 
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where 



(B18) 



vl = hl + hl- 2hih^a2,i, (B19) 



«2 



«:^ = l-«2^-(l + 4J^, (B20) 

^0 

a+ = -(hl-hl)^. (B21) 
Combining Eqs. ( |Bl(j| ) and ( [BIT] ) yields the solutions 

"3°^ (r) = as^i cos $ + //(a+^j - a+) sin $ (B22) 
a+(r) = d+ sin $ + — a+) cos $, (B23) 

where a+,i = a+l'"*), "s.i = a3(n), and 

$ = z/o(r-ri). (B24) 

Note that uq is the frequency of oscillation of the af^^s — frequency with respect to changing 
orbital radius r, not changing time t. 

Analytic expressions for the functions ai^\ . . . , ai^"* can now be obtained by combining 
Eqs. ( [B13D - ( |B15| ) and ( |B18| ) - ( |B23| ). These expressions depend in a complex way on all 
of the initial spin direction parameters ai^i, a^^i and 04^4, and also on S2, on the mass ratio 
M1/M2, and on the initial orbital separation r^/M. For the equal mass case Mi = M2, the 
"frequency" z/q is given by 



s\ r 255 1125 



On 



4096 ■ 16384'^^'*. 



(B25) 



Values of for Mi 7^ M2 are typically much larger than this. 

The first-order corrections OL'^^ir') can be obtained using the zeroth order solutions and 
Eqs. ( p6|) - ( |B8|) . However, we are only interested in determining the leading order behavior 



of /3, and for this purpose we need only evaluate 0^2^ . From Eqs. ( P7| ) and ( |B22| ), this is 
given by 
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15 



128;Uz/o 



Wi 



L 



X [0:3 j sin $ — z/ (a 



ce^) (cos<l> - 1)1 . 



(B26) 



Substituting Eqs. 



and ( |B18| ) - ( |B23| ) into ( |B9D gives a result of the form 



/5(r) = A + i? cos $ + C sin <l>, 



(B27) 



where the constants B and C are first order in Si. The resulting expression for the total 
amplitude of oscillation A = ^/ B^ + in terms of the variables ai^i, . . . , a^^i, M1/M2 and 
r is complicated and not very illuminating, so we do not reproduce it here. Instead we show 
in Figs. and ^ the quantity ^max = -^max[«2,j, 04,1] obtained in the following way: (i) 
Use Eq. ( [B5| ) to eliminate a^^i in terms of ai^i, 02,4 and a^^i. (ii) Numerically maximize over 



values of ai.j that lie in the range between the values a2,i ce^^i ± y (1 — a| J (1 — a;| J allowed 
by Eq. (B3). (iii) Choose the maximal spin magnitudes Si = Mf, S2 = M|. (iv) Choose 
the final orbital separation r = 6M, the value for which the amplitude A will most likely be 
largest. It can be seen from these plots that for all choices of initial angles, A < 0.25. 

In the special case that Mi = M2, the formulae simplify and we find that A oc 
(this is not true in general). Specifically we find in this case that 

376 



B 



OiiA a2A - 2a4 , 



384v/F 



(B28) 



and 



C 



235 a2,i«3,; 



(B29) 



512 v^z/o ' 

where Vq is given by Eq. ( B25 ). 

To summarize, we have determined the evolution of the quantity x = 47r — /3 both 
numerically, for a wide range of initial conditions, and analytically, in the regime where 
|S^i| ^ M^. In all cases we find that the amplitude of the oscillations of x is < 0.25. 
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APPENDIX C: THE DECOUPLING OF PHASE AND 
AMPLITUDE PARAMETERS IN THE FISHER 
INFORMATION MATRIX 



The phase ^'(/) of the Fourier transform of the waveform can be written in the form 



n.=l,2,--- 



(CI) 



where (ai, 02, as, • • •) = (O5 ^5/3, —1, —2/3, . . .), and the parameters ci, C2, C3, C4 etc. are 
simply related to the parameters 0c, ^c, A^, /^, etc., via Eq. (|5BD. The number of variables 
Cn will depend on the post-Newtonian order to which \E'(/) is calculated; the following 
analysis holds for any number of these variables. We can make a linear transformation to 
new variables 



d =11 

"1 m 



(C2) 



in such a way that 



and that for m > 2 |^ 
dh 



dh dh 



ih 



ddi dci 



roo l/)nl^ 



ih. 



(C3) 



m/ ri 

0. 



(/) 



(C4) 



The key point now is that the inner product 



dh 

dfi" 



dh 

ddm. 



where fi"' is any of the "amplitude" parameters D, ip, v and di, will also be proportional to 
the right-hand side of Eq. (|C^ ) for m > 2, and so will vanish. This can be seen from the 
structure of Eq. (]66|). Consequently, in the new variables D, v, ip, and dm, m = 1,2, . . ., the 
Fisher matrix ( [A31| ) will be block diagonal, which establishes the result stated in Sec. [IV C 
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FIGURES 

FIG. 1. Gravitational waveforms from coalescing compact binaries are completely specified 
by a finite number of parameters 9 = {9^, . . . ,9^), and so form a surface S in the vector space 
V of all possible measured detector outputs s = s{t). The statistical properties of the detector 
noise endow V with the structure of a infinite-dimensional Euclidean space. This figure illustrates 
the relationships between the true gravitational wave signal h{9), the measured signal s, and the 
"best-fit" signal h{9). Given a measured detector output s = h{9) + n, where n = n{t) is the 
detector noise, the most likely values 9 of the binaries parameters are just those that correspond 
to the point h(9) on the surface S which is closest [in the Euclidean distance {s — h \ s — h)] to y. 

FIG. 2. This plot shows how the total signal-to- noise squared S'^ /N'^ for a detected coalesc- 
ing-binary waveform is distributed in frequency /, assuming the detector noise curve @. Most of 
the signal-to-noise ratio comes not near 70 Hz where the detector sensitivity Sn{f)~^ is highest, 
but rather at a somewhat lower frequency of ~ 50 Hz, because more cycles per unit frequency are 
received at lower frequencies. 

FIG. 3. This plot shows the curve of constant probability on the fi(3 plane for a NS-BH binary, 
where fi is the binary's reduced mass and /3 is a dimensionsless spin-related parameter, such that the 
true values of these parameters lie inside the ellipse with 95% confidence. The strong correlation 
between possible values of fx and /3 is evident. To a good approximation, the chirp mass M is 
measured to arbitrarily high accuracy. Hence in the three-dimensional space of the parameters 
{M,fj,,l3), the true values of these parameters are confined with high confidence to a thin strip of 
the above ellipsoidal shape in a plane of constant Ai. 
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FIG. 4. A diagram showing the information obtained from the gravitational wave signal, con- 
straining the individual masses Mi and M2 of the binary components, in various cases. Because 
of the highly accurate measurement of the chirp mass M in each case, the individual masses are 
essentially constrained to lie on a curve of constant M = {MiM2f'^{Mi + Ms)-^/^ in the Mi M2 
plane. The measured value of Ai provides a sharp lower bound for the larger mass Mi, and a 
sharp upper bound for the smaller mass M2, since the constant- curves terminate sharply at the 
forbidden, hatched region. The measurement of the reduced mass /x gives some (but not much) 
information about where along the constant Ai curve the binary is most likely to be located. In 
each case, the solid circles show the true values of Mi and M2, the solid curve denotes the 68% (1 
sigma) confidence interval, and the dashed extension denotes the 95% (2 sigma) confidence interval. 
The detector noise spectrum (Q) was assumed. 

FIG. 5. The amplitude sensitivity function cj£)(n), as a function of position on the sky 
parametrized by the Earth-fixed coordinates 9 and ip, for the detector network consisting of the 
two LIGO detectors in Hanford, Washington and Livingston, Louisiana, and the VIRGO detector 
in Pisa, Italy. The axis ^ = is the Earth's axis of rotation, and (/? = is 0° longitude. Only sky 
positions over the northern hemisphere are shown, because an takes the same values at antipodal 
points. The function cr/3(n) has the following meaning: for a source of waves in the direction n, 
the combined signal-to-noise ratio of the whole network, averaged over all polarization angles tp of 
the source (equivalently, averaged over rotations of the source in the plane perpendicular to the 
line of sight), will be proportional to aoin). The thick black lines indicate the positions of the 
three detectors. This plot can be generated by combining Eqs. ( |8^ ) and ( pO[ ) of the text with the 
network parameters given after Eq. ([90|). 
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FIG. 6. The polarization sensitivity function 1— 5£)(n), for the LIGO/VIRGO detector network; 
see caption of Fig. |5|. This plot can be generated by combining Eqs. (|85[), (^), (|89|), and ( |90| ) of the 
text. The quantity ez)(n) essentially measures the "skewness" or assymetry in the sensitivities of 
the network to the two independent polarization components of waves propagating in the direction 
n. When ~ 0, the network has roughly equal sensitivity to both polarization components. When 
£d ~ 1, on the other hand, one polarization component can be measured far more accurately than 
its orthogonal counterpart. In this case the signal-to-noise ratio for incident, strongly linearly 
polarized bursts of waves (e.g., those from edge-on coalescing binaries) will depend sensitively on 
the polarization axis, i.e., it would vary by large factors if the source were rotated in the plane 
perpendicular to the line of sight. Note that the polarization sensitivity is poor 0.2) for directions 
directly overhead the two LIGO detectors (because the two detectors are nearly parallel), and is 
typically ^ 0.3 over most of the sky. Good sensitivity is achieved in isolated regions. 

FIG. 7. The dependence of the distance measurement accuracy /S.D/D on the sky location n, 
the polarization angle if), and the cosine v of the angle of inclination of the orbit to the line of sight 
is approximately given by /S.D/D oc T(n, t;,^/;), where the dimensionless function T is defined in 
Eq. (|97|). Here we plot for the LIGO/VIRGO detector network the quantity Tmax obtained by 
maximizing T over all polarization angles if), at = 1/2, as a function of 6 and (p. Higher values 
of T indicated by regions of lighter shading correspond to poorer measurement accuracy. 

FIG. 8. The quantity Tmin which is obtained by minimizing T(n, v,-)/^) over at v'^ = 1/2; see 
caption of Fig. |^ 



96 



FIG. 9. The quantity i}{eD), which is the sohd angle on the sky for which the polarization 
sensitivity is less than 1 — en, for two different detector networks. The solid line indicates the 
LIGO/VIRGO detector network, and the dashed line a 4-detector network consisting of the LIGO 
and VIRGO detectors together with a hypothetical detector in Perth, Australia. These plots 
were generated using 1000 randomly chosen directions n. The great improvement in polarization 
sensitivity due to the additional detector is apparent: e.g., the polarization sensitivity is < 0.2 over 
~ 60% of the sky for the 3-detector network, but only over ~ 20% of the sky for the 4-detector 
network. 

FIG. 10. An example illustrating the necessity of going beyond the Gaussian approximation. 
Consider a neutron-star neutron-star binary merger in the direction given by {6,<f) = (50°, 276°). 
The LIGO/VIRGO network parameters for this direction are an = 1-03 and ed = 0.74. Sup- 
pose that an experimenter determines from the measured signal the following "best-fit" (maxi- 
mum-likelihood) parameters: distance Dq = 432 Mpc [corresponding to a signal-to-noise ratio of 
p = 12.8, assuming the advanced detector sensitivity level (^)], masses Mi = M2 = 1.4Mq, cosine 
of inclination angle vq = 0.31, and polarization angle ipo = 56.5°. Then the distribution that she 
would infer by a Bayesian analysis for the distance to the source is shown by the solid curve; it 
is given by Eqs. (|112| ), (|115| ), and (^) of the text. The Gaussian approximation [Eq. (|96|) of the 
text] to this distribution is shown by the dashed curve. The distance measurement accuracy is 



atypically poor in this example; see Fig. 16 below. 



FIG. 11. The solid line shows the distance measurement accuracy AD/D for the binary merger 
discussed in the caption of Fig. |l^ (for which Ai = 0.10, A2 = 0.057), as a function of the cosine 
of the angle of inclination, vq. The improvement in accuracy at high values of vq is due in part to 
an increased signal-to-noise ratio there. The dashed curve shows the prediction ( |96[ ) of the linear 
error-estimation theory, which diverges as fo — > 1- 
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FIG. 12. The distance measurement accuracies that result from displacing along the line of 



sight to the Earth, to various distances Dq, the binary merger of Fig. |10|. As in Fig. 11, the dashed 
curve shows the approximate hnear estimate (|9^), and the solid curve shows the more accurate 
estimate ( [1171) . The curves terminate at that distance (~ 700 Mpc) beyond which the merger is no 
longer visible, assuming the detector sensitivity level (|^ and a combined signal-to-noise threshold 
of 8.5. 

FIG. 13. The distances Dq for 1000 NS-NS binaries whose locations and orientations were 
randomly chosen, and the corresponding predicted signal-to-noise ratios p for the LIGO/VIRGO 
network. The lower dashed line is the signal-to-noise threshold of 8.5, below which sample points 



were discarded; the upper dashed line shows the maximum possible value (104) of p at a given 
distance. Six points with p between 50 and 90 are not shown. The number of sources with p larger 
than a given value p^: is proportional to p^'^. The detection of this many binary inspirals with 
the advanced LIGO/VIRGO detectors would take several years, if merger rates are as currently 
estimated MM , and assuming the detector sensitivity level (0) . 



FIG. 14. The distance measurement accuracy AD/D computed from Eq. ( |117| ) for the same 
1000 NS-NS binaries, versus the distance Dq; see caption of Fig. |13[ The spread in the values of 
AD/D is due to different source directions and orientations. Note that the accuracy for sources 
within 200 Mpc (of which there are estimated to be ~ 3 per year Ij^Q) can vary between ~ 2% 
and ~ 25%. For the most distant detectable sources (at ~ 1200 Mpc), the accuracy can sometimes 
be as good as ~ 20%. The dashed line shows the theoretical lower bound ( |102D derived using the 
linear error-estimation formalism; points below this line mostly have values of vq close to one for 
which value the linear error-estimation theory fails. 

FIG. 15. Distance measurement accuracy versus signal-to-noise for the same 1000 NS-NS bi- 
naries; see caption of Fig. The dashed line shows the theoretical lower bound ( |103| ). As in 
Fig. O, six points with p between 50 and 90 are not shown. 
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FIG. 16. The frequency of occurrence of different ranges of AI? /D, out of a total of 1000 signals, 
for the LIGO/VIRGO detector network. It can be seen that ~ 8% of detected signals will have 
distance measurement accuracies of better than 15%, while ~ 60% of them will have accuracies of 
better than 30%. These conclusions are insensitive to the overall scale of the detectors' intrinsic 
noise, which essentially sets the event-detection rate. By contrast, they are sensitive to the number 
of detectors in the detector network, and to their orientations; see text and also Fig. ^. 

FIG. 17. Distance measurement accuracy versus the detector network polarization sensitivity 



1 — e£)(n), for 1000 NS-NS binaries; see caption of Fig. 13. The strong correlation between very 
poor distance-measurement accuracy and low polarization sensitivity is evident — essentially all 
points with AD/D > 0.5 have 1 - < 0.2. 

FIG. 18. During the last few minutes of inspiral, the angles between the bodies' spins ^i, S2 
and the orbital angular momentum L all oscillate, in addition to and separately from the precession 
of the total spin S = Si+ S2 around L. This oscillation gives rise to an oscillation of the parameter 
X = 47r — /3 which governs the contribution of the spins to the accumulated phase of the emitted 
gravitational waves [cf. Eq. (|49|)]. Here we show a typical plot of x as a function of the orbital 
separation r, for a NS-NS binary. The spin and orbital angular momentum directions were taken 
to be 5*1 oc i + k, S2 oc —j and L oc i + j at the initial gravitational wave frequency of 10 Hz. The 
spin of each neutron star was assumed have the maximal magnitude of (1.4Mq)^, corresponding to 
a rotation period of a few miliseconds. It can be seen that in each case the amplitude of oscillation 
of X is very small compared to its mean value, so that to a good approximation we can take 
X = constant. 



FIG. 19. The evolution of x ^or a NS-BH binary; see caption of Fig. 18. The black hole was 
assumed to be maximally rotating. 



FIG. 20. The evolution of x for a BH-BH binary; see caption of Fig. 18. Both black holes were 
assumed to be maximally rotating. 
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FIG. 21. When one of the spins is small, the evolution of the parameter x = 47r — /3 is ap- 
proximately given by X = Xo + cos [vor + const] , where the "frequency" uq and amplitude A 
are slowly varying functions of r. Here we show the amplitude A for an equal-mass binary, at an 
orbital separation of r = 6M, as a function of 02, i = 5*2 • and j = Si - S2, where the maximum 
is taken over the remaining angles. Si, S2 and L are unit vectors in the directions of the initial 
spins and the initial orbital angular momentum. 



FIG. 22. As in Fig. 21, but for a NS-BH binary with M1/M2 = 1.4/10. 
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TABLES 

TABLE L The rms errors for signal parameters and the correlation coefficient CMfi, calculated 



assuming spins are neg 


ligible. The results are for 


a single "advanced" detector, the shape of whose 


noise curve is given by Eq. 


Ml and M2 are 


in units of solar masses, 


while At/; is in 


units of 


msec. 


The rms errors 


are normalized to a signal-to-noise ratio of S/N = 


= 10; the errors 


scale as 


(S/N) 


while CMfi is 


independent of S/N. 








Ml 


M2 




At, 


AM/M 


All / II 




2.0 


1.0 


1.31 


0.721 


0.0038 % 


0.39 % 


0.899 


L4 


1.4 


1.28 


0.713 


0.0040 % 


0.41 % 


0.906 


10 


1.4 


1.63 


1.01 


0.020 % 


0.54 % 


0.927 


15 


5.0 


2.02 


1.44 


0.113 % 


1.5 % 


0.954 


10 


10 


1.98 


1.43 


0.16 % 


1.9 % 


0.958 
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TABLE II. The rms errors for signal parameters and the correlation coefficients cjUfi, cm/b, and 
c^^, calculated using spin-dependent waveforms. The results are for a single "advanced" detector, 
the shape of whose noise curve is given by Eq. (^). For the rows marked with a f (and only for 
those rows), the variance-covariance matrix has been "corrected" to approximately account for the 
fact that the spin parameter /3 must satisfy < (3max ~ 8.5. The rms errors are normalized to a 
signal-to-soise ratio of S/N = 10. Except for rows marked with a f, errors scale as (S/N)^^, while 
the correlation coefficients are independent of S/N. Except for rows marked with a f, if /? had 
been chosen non-zero with Mi and M2 unchanged, then AAi/Ai, A/i/7W, and cm^i would have 
been unchanged (but A/?, c^k/j, and c^/j would have been altered). As in Table |, Mi and M2 are 
in units of Mq, while Ate is in msec. The results for the LIGO /VIRGO network of detectors, for 
a signal with combined signal-to-noise ratio from all the detectors of 10, will be approximately the 



same 


as those shown here; 


see text. 
















Ml 


M2 


13 


A^e 


Ate 


AM/M 


Afi/ij. 


A/3 


CM^l 


CMP 




2.0 


1.0 





4.13 


1.14 


0.034 % 


8.44 


% 


1.04 


-0.988 


0.993 


-0.9989 


1.4 


1.4 





4.09 


1.13 


0.034 % 


9.65 


% 


1.24 


-0.988 


0.993 


-0.9991 


10 


1.4 





6.24 


2.03 


0.19 % 


15.2 


% 


1.99 


-0.990 


0.994 


-0.9994 


5 


1.4 





4.89 


1.44 


0.10 % 


13.4 


% 


1.73 


-0.989 


0.994 


-0.9992 


15 


5 





9.26 


3.53 


1.06 % 


76.4 


% 


11.4 


-0.992 


0.994 


-0.99980 


15 


5 


ot 


5.77 


2.40 


0.64 % 


45.8 


% 


6.81 


-0.978 


0.984 


-0.9995 


10 


10 





9.26 


3.53 


1.42 % 


125 \ 




19.5 


-0.992 


0.994 


-0.99988 


10 


10 


ot 


4.13 


1.92 


0.59 % 


49.9 


% 


7.79 


-0.953 


0.964 


-0.9992 
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TABLE in. Measurement accuracies, including spins, as in Table || except that we the take 



shape of the 


noise 


curve to b 


e given 


by Eq. (g. 












Ml 


M2 


(3 




Ate 


AM/M 


Afi/fi 


A/3 




CMI3 


Cm/3 


2.0 


1.0 





2.21 


0.47 


0.021 % 


5.12 % 


0.641 


-0.986 


0.992 


-0.9988 


1.4 


1.4 





2.19 


0.46 


0.021 % 


5.84 % 


0.757 


-0.986 


0.992 


-0.9990 


10 


1.4 





3.90 


1.06 


0.13 % 


10.1 % 


1.33 


-0.990 


0.994 


-0.9993 


5 


1.4 





2.81 


0.66 


0.065 % 


8.47 % 


1.11 


-0.988 


0.993 


-0.9992 


15 


5 





6.72 


2.40 


0.75 % 


55.1 % 


8.18 


-0.992 


0.994 


-0.99980 


15 


5 


Of 


4.95 


1.85 


0.55 % 


39.7 % 


5.90 


-0.985 


0.989 


-0.9996 


10 


10 





6.72 


2.40 


1.00 % 


90.0 % 


14.0 


-0.992 


0.994 


-0.99988 


10 


10 


ot 


3.71 


1.49 


0.53 % 


46.7 % 


7.27 


-0.972 


0.978 


-0.9995 
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